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Time  series  models  with  autoregressive,  moving  average 
and  mixed  autoregressive-moving  average  correlation  struc¬ 
ture  and  with  positive-valued  non-normal  marginal  distribu¬ 
tions  are  considered.  First, a  flexible  mixed  model 
GLARMA (p , q)  with  Gamma  marginals  is  investigated.  The 
correlation  structure  for  several  special  cases  is  derived. 
For  the  first-order  autoregressive  case,  GLAR(l),the 
conditional  density  of  X  given. X  .  is  derived.  This 
leads  to  the  formation  of  a  likelihood  function  and  a 
numerical  approximation  to  and  a  simulation  study  of  the 
maximum  likelihood  method  of  parameter  estimation.  Multi¬ 
variate  extensions  of  the  model  are  considered  briefly. 

Second,  three  methods  for  generating  first-order 
moving  average  sequences  with  Exponential  marginals  are 
examined.  These  generalize  the  EMA(l)  Exponential  model. 
Negative  correlation  using  antithetic  variables  is  investi¬ 
gated  in  the  moving  average  models. 

A  preliminary  analysis  of  wind  speed  data  obtained  over 
a  15  year  period  in  the  Gulf  of  Alaska  is  presented.  A 
model  with  four  harmonic  deterministic  mean  multiplying 
random  innovative  factors  modeled  by  a  GLAR(l)  process  is 
developed.  Correlograms  and  periodograms  are  used  to  deter¬ 
mine  the  model  for  the  mean  and  the  structure  of  the 
innovation  process. 
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IV.D.6g  Correlogram  1961  detrended  data  -  314 

IV.D.6h  Correlogram  1962  detrended  data  -  315 

IV.D.6i  Correlogram  1963  detrended  data  -  316 

IV. D.  6  j  Correlogram  1964  detrended  data -  317 
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IV.D.6k 
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IV. D 
IV. D 
IV. D 
IV. D 
IV. D 
IV. D 
IV. E 

IV. E 

IV. E 

IV. E 

IV. E 

IV. E 

IV. E 

IV. F 

IV.  F 

IV.  F 

IV. G 

IV. G 


Correlogram  1965  detrended  data 

61  Correlogram  1966  detrended  data -  319 

6m  Correlogram  1967  detrended  data -  320 

6n  Correlogram  19  68  detrended  data -  321 

6o  Correlogram  1969  detrended  data -  322 

6p  Average  correlogram,  detrended  data  -  323 

7  Correlogram  of  average  detrended  data  -  324 

1  Average  data  plotted  against  2 

harmonic  smoothed  data  -  326 

2  Periodogram  1955  2  harmonic  detrended 

data  with  AR  ( 1)  spectrum -  327 


3  Log  periodogram  1955  2  harmonic 

detrended  data  with  log  AR(1)  spectrum  —  328 

4  Periodogram  1969  2  harmonic  detrended 


data  with  AR(1)  spectrum -  329 

5  Log  periodogram  1969  2  harmonic 

detrended  data  with  log  AR(1)  spectrum  —  330 

6  Periodogram  15  year  average  2 

harmonic  detrended  data  -  333 

7  Log  periodogram  15  year  average 

2  harmonic  detrended  data  -  334 

1  Periodogram  prewhitened  15  year 

average  2  harmonic  detrended  data  -  336 

2  Log  periodogram  prewhitened  15  year 

average  2  harmonic  detrended  data  -  337 

3  Correlogram  prewhitened  15  year  average 

2  harmonic  detrended  data  -  338 

1  Periodogram  1955  4  harmonic  detrended 

data  with  AR(1)  spectrum -  341 

2  Log  periodogram  1955  4  harmonic 

detrended  data  with  log  AR(1)  spectrum  —  342 


15 


IV. G. 3 


Periodogram  1969  4  harmonic  detrended 
data  with  AR (1)  spectrum  - 


343 


IV. G. 4 

IV. G. 5 

IV. G. 6 

IV.  G.  7 


Log  periodogram  1969  4  harmonic 
detrended  data  with  log  AR(1)  spectrum  —  344 

Periodogram  15  year  average  4  harmonic 
detrended  data -  345 

Log  periodogram  15  year  average 
4  harmonic  detrended  data  -  346 

GLAR(l)  Sample  path;  k=4.0,  q=0.6;  1 

TRUE  Rho:  0.85  -  348 


1 
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A 


I .  INTRODUCTION 


The  classical  approach  to  time  series  analysis  based  on 
linear,  additive  models  with  normally  distributed,  constant 
variance  residuals  is  probably  best  presented  by  the  work  of 
Box  and  Jenkins  [Ref.  1] .  Although  their  work  is  widely  ac¬ 
cepted  and  used,  it  is  not  applicable  to  some  important  time 
series.  This  is  mainly  because  the  Box-Jenkins  approach  is 
based  on  an  assumed  normal  distribution  for  the  series  in 
question.  However,  the  assumption  of  normality  is  not  appro¬ 
priate  when  the  series  is  known  to  be  non-negative.  Such 
series  typically  involve  times  between  successive  events  in 
event  processes.  Examples  are  easy  to  construct.  Times  be¬ 
tween  arrivals  at  a  hospital  emergency  room,  times  between 
breakdowns  in  a  tank  main  drive  assembly,  and  times  between 
detections  of  enemy  armor  vehicles  are  a  sample  of  series  of 
this  type.  Because  of  the  non-negative  nature  of  the  series, 
the  Box-Jenkins  distributional  assumptions  and,  hence,  the 
analysis  techniques  are  inappropriate.  There  is,  of  course, 
the  possibility  of  data  transformations  but  this  is  not  appro¬ 
priate  with  very  skewed  marginal  distributions  and  it  is,  in 
most  cases,  difficult  to  ascertain  what  the  transformation 
does  to  the  correlation  structure  of  the  series. 

Gaver  and  Lewis  [Ref.  2]  wrote  the  pioneering  paper  on 
the  subject  of  autoregressive  processes  with  non-normal  marginal 
distributions.  They  presented  the  method  for  determining  the 
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distribution  of  the  innovative  terms  in  the  basic,  linear, 
additive,  autoregressive  equations  (first-order  stochastic 
different  equation) 


X 


n 


PX  ,  +  £ 
n-1  n 


(1.1) 


that  was  required  to  produce  a  given  marginal  distribution 
for  the  {Xn>  sequence.  They  presented  results  for  {Xfl}  se¬ 
quences  with  Exponential,  Gamma,  and  mixed  Exponential 
marginals.  They  also  showed  that  this  problem  was  the  same 
as  that  of  determining  the  class  of  self-decomposable  (Type 
L)  random  variables  (Feller,  [Ref.  3],  Loeve  [Ref.  4])  al¬ 
though  the  connection  between  the  solution  to  the  self-decom- 
posable  problem  and  equation  (1.1)  was  not  explicit  in  the 
literature. 

The  Gaver  and  Lewis  paper  was  followed  by  other  papers 
which  extended  these  results.  Lawrance  and  Lewis  [Ref.  5] 
presented  a  first-order  moving  average  process  with  Exponential 
marginals.  Jacobs  and  Lewis  [Ref.  6]  propounded  a  mixed  auto¬ 
regressive-moving  average  of  order  one,  EARMA (1,1) ,  and 
Exponential  marginals.  This  was  extended  to  an  arbitrary 
order  EARMA(p,q)  process  by  Lawrance  and  Lewis  [Ref.  7],  A 
further  refinement  of  the  first-order.  Exponential,  auto¬ 
regressive  process  (NEAR(l) )  was  presented  by  Lawrance  and 
Lewis  [Ref.  8] .  While  this  contained  the  previous  EAR(l) 
model,  it  did  not  suffer  from  the  degeneracy  inherent  in  (1.1). 
Jacbos  applied  these  models  to  closed  cyclic  queueing  networks 
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[Ref.  9]  and  Lewis  and  Shedler  applied  them  to  models  of 
computer  processes  [Ref.  10]. 

This  paper  extends  the  results  of  these  researchers  and 
others  in  three  areas.  In  Chapter  II  a  mixed  ARMA(p,q)  model 
with  Gamma  marginals  proposed  by  Lewis  [Ref.  11] ,  the  GLARMA(p,q) 
model,  is  examined.  The  correlation  structure  is  derived  for 
several  values  of  p  and  q.  Of  particular  note  is  the  AR(1) 
case  (p  =  1,  q  =  0)  ,  called  GLAR(l) ,  where  the  conditional 
density  of  Xn  given  Xn_-^  is  derived.  This  leads  to  the  deri¬ 
vation  of  a  likelihood  function  and  a  numerical  technique  to 
evaluate  and  maximize  the  likelihood  function  with  respect  to 
the  model  parameters.  This  provides  a  useful  technique  for 
estimating  model  parameters.  Using  this  numerical  technique 
a  simulation  study  of  the  properties  of  maximum  likelihood 
estimators  for  the  parameters  of  the  model  is  given. 

The  correlation  structure  is  derived  for  other  models  in 
the  GLARMA ( p , q )  family:  the  first-order  moving  average,  the 
second-order  autoregressive,  the  first-order  mixed  autoregres¬ 
sive-moving  average  and  a  bivariate  first-order  autoregressive 
process.  These  different  models,  particularly  the  bivariate 
extension,  demonstrate  the  flexibility  of  the  GLARMA (p,q) 
model. 

In  Chapter  III  the  first-order  moving  average  process  with 
Exponential  marginals  of  Lawrance  and  Lewis  [Ref.  5]  is  ex¬ 
tended  to  a  two  parameter  model.  This  is  done  by  utilizing 

\ 

the  NEAR { 1 )  structure  which  combines  two  independent  Exponen¬ 
tial  random  variables  into  a  random  variable  with  Exponential 

< 
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distribution.  A  fairly  complete  set  of  characteristics  of 
this  model  are  derived.  In  particular  the  correlation  struc¬ 
ture,  the  quantity  P(Xn+^>Xn),  the  Laplace  transform  of 
sums  of  Xn's,  the  Laplace  transforms  of  the  distribution  of 
counts,  the  (Bartlett)  spectrum  of  counts,  and  the  joint 
Laplace-Stielt jes  transforms  of  Xn  and  Xn+^  are  addressed. 

These  characteristics  are  compared  to  those  of  other  proc¬ 
esses  which  produce  marginally  Exponential  random  variables. 

In  Chapter  IV  the  models  of  Chapter  II  are  used  in  a 
preliminary  data  analysis  of  wind  speed  data.  This  repre¬ 
sents  the  first  effort  to  apply  these  models  to  a  large,  real 
world  data  base.  A  model  for  simulating  wind  data  is  pre¬ 
sented  and  parameter  estimates  for  the  data  are  derived  using 
the  numerical  approximation  to  the  maximum  likelihood  presented 
in  Chapter  II. 


II.  MODELS  WITH  GAMMA  MARGINALS 


A.  INTRODUCTION 

There  have  been  several  scheir.es  suggested  for  the  modeling 
of  dependent  random  variables  with  Gamma  marginals.  The 
Gamma  autoregressive  process  of  order  one  (GAR(l) )  by  Gaver 
and  Lewis  [Ref.  2],  the  discrete  autoregressive  process  of 
order  one  (GDAR(l))  by  Jacobs  and  Lewis  [Ref.  12],  the  Gamma 
Beta  autoregressive  process  of  order  one  { GBAR ( 1 ) )  by  Fishman 
[Ref,  13]  and  Lawrance  and  Lewis  [Ref.  14],  and  the  Gamma  auto¬ 
regressive  process  of  order  one  ( GLAR  ( 1 )  )  by  Lewis  [Ref.  10]. 
There  is  also  an  attempt  to  use  multivariate  Gammas  obtained 
by  the  inverse  probability  integral  transform  in  a  time  series 
context  by  chmeiser  [Ref.  15]  . 

The  GAR(l)  model  generates  an  {X^}  series  using  the  stan¬ 
dard  first-order  autoregression  equation  (first-order  stochas¬ 
tic  difference  equation) 

Xn  "  oXn-l  +  £n'  0  1  0  <  (II. A. 1) 

The  innovative  factor,  en,  has  Laplace-Stielt jes  transform 

\  k 

of  (p+(l-p)-r^)  and  the  i. Xn }  random  variables  have  marginal 
Gamma  distributions  with  shape  parameter  k  and  scale  parameter 
A.  The  marginal  density  function  of  the  {Xn>  random  variable 


fx  (x;k,  \)  =  fx(x;k,X)  =  r^y-  xk_1  e~Ax,  (II. A. 2) 

A.>0,x>0,k>0. 

The  model  tends  to  produce  runs  of  decreasing  value  when  inno¬ 
vative  term  has  successive  realizations  of  value  zero.  The 
GAR ( 1)  is  in  this  sense  highly  degenerate,  even  though  it  is 
a  true  linear  process.  Ad  hoc  estimates  of  model  parameters 
are  available  which  produce  the  exact  p  value  if  the  series 
is  long  enough  [Ref.  2] .  However,  maximum  likelihood  esti¬ 
mates  have  not  been  produced.  This  model  is  not  extendable 
to  a  moving  average  process. 

The  GDAR(l)  produces  an  { }  sequence  using  the  first- 
order  autoregressive  equation  with  random  coefficients. 


n 


=  V  X  .  +  (1-V  )G  , 
n  n- 1  n  n 


(II. A. 3) 


where  {Vn,  n  =  1,2,...}  is  an  iid  sequence  of  binary  random 

variables  with  P(V  =1)  =  1-?(V  =0)  =  p,  (G  ,  n  =  1,2,...; 

n  n  n 

is  an  iid  Gamma  sequence. 

This  sequence  produces  runs  of  constant  value  when  suc¬ 
cessive  realizations  for  V  produce  value  1.  When  V  equals 

n  n 

zero,  a  new  value  is  selected.  Obviously,  this  model  has 
limited  value  in  general  applications  and  is  even  more  degener¬ 
ate  than  GAR ( 1)  process. 

The  GBAR(l)  is  the  most  flexible  model  in  that  it  contains 
the  GAR(l)  and  GLAR(l)  models  as  special  cases.  It  produces 


24 


an  {X^}  sequence  using 


xn  "  6Bnxn+l  +  S'  (II-A"» 

where  {Bn,  n  =  1,2,...}  is  an  iid  Beta  (k-q,q)  sequence.  en 
was  shown  by  Lawrance  and  Lewis  to  be  the  sum  of  a  Gamma 
variable  and  the  innovation  process  of  the  GAH(l)  process 
[Ref.  14]  .  Although  flexible  in  the  sense  that  it  contains  the 
other  models,  it  can  not  be  extended  to  a  moving  average  proc¬ 
ess.  In  addition,  conditional  densities  and,  hence,  maximum 
likelihood  estimates  are  not  available.  This  is  because  the 
innovation  random  variable  for  the  GAR(l)  process,  while  it 
can  be  generated  as  a  random  sum  of  random  variables,  does  not 
have  a  known  distribution  function. 

The  most  valuable  and  flexible  model  seems  to  be  the 
GLAR(l)  which  produces  an  {X^}  sequence  using  the  stochastic 
difference  equation  with  random  coefficients 


X. 


n 


=  B  X  .  +  C  G  , 
n  n-1  n  n 


(II. A. 5) 


where  (Xn,  n  =  0,1,...}  is  a  second-order  stationary  sequence 
of  Gamma  random  variables,  {B  ,  n  =  1,2,...}  and  { ,  n  =  1,2,... 
are  iid  Beta  random  variables,  and  {Gn}  is  an  iid  sequence  of 
Gamma  random  variables.  This  model  is  extendable  as  an  auto¬ 
regressive  process  of  arbitrary  order  and  as  a  moving  average 
process  of  arbitrary  order.  These  two  forms  can  also  be  combined 
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to  form  a  mixed  model,  the  so-called  Gamma  Lewis  autoregressive  - 
moving  average  process  of  orders  p  and  q  (GLARMA (p, q)  )  . 

This  chapter  of  the  thesis  examines  some  of  the  special 
cases  of  the  model.  One  case  in  particular,  the  AR(1)  form, 
is  reasonably  extensively  examined.  The  correlation  struc¬ 
ture  is  developed.  The  conditional  expectation  and  density 
are  derived.  The  latter  is  used  as  the  basis  of  a  numerical 
approximation  to  the  maximum  likelihood  method  of  parameter 
estimation.  Directional  moments  and  the  probability  of  Xn+1 
being  greater  than  X^  are  derived.  In  a  later  Chapter  of 
this  thesis,  this  model  is  used  as  a  basis  for  analyzing  wind 
speed  data. 

The  special  case  of  the  moving  average  of  order  one  is 

examined  in  some  detail.  The  correlation  structure  is  derived 

with  some  emphasis  on  exploring  the  restrictions  on  the  range 

of  correlations  that  are  possible.  Directional  moments  and 

an  empirical  examination  of  the  probability  that  Xn+^  is 

greater  than  X  are  examined, 
n 

As  a  demonstration  of  the  flexibility  and  extendabilitv 
of  the  model  the  mixed  model  of  order  one,  the  autoregressive 
model  of  order  two,  and  a  bivariate  model  are  introduced  and 
their  correlation  structures  derived. 

B.  FIRST-ORDER  AUTOREGRESSIVE  BETA-GAMMA  MODEL,  GLAR(l) 

1 .  Introduction 

The  first-order  autoregressive  Beta-Gamma  model  is  a 
special  case  of  the  GLARMA (p,q)  model  when  q  =  0  and  p  =  1. 
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The  autoregressive  model  generates  an  {Xn}  sequence  using 

xn  ’  Anx„-1  +  BnV  HI. a. 1.1) 

where  {Xn,  n  =  0,1,...}  is  a  second-order  stationary  sequence 
of  random  variables  with  a  Gamma (k,l)  marginal  distribution; 

{ An ,  n  =  1,2,...}  is  an  iid  sequence  of  Beta(k-q,q)  random 
variables;  {B  ,  n  =  1,2,...}  is  an  iid  sequence  of  Beta (q,k-q) 
random  variables;  {Gn,  n  =  0,1,...}  is  an  iid  sequence  of 
Gamma  (k,l)  random  variables;  {An} ,  {Bn},  and  {Gn}  are  inde¬ 
pendent;  0  <  q  <_  k. 

Choosing  XQ  =  GQ  makes  the  {Xn}  sequence  stationary. 

The  Gamma  random  variables  are  parameterized  by  the  shape 
parameter  and  the  mean,  rather  than  the  scale  parameter.  This 
somewhat  unusual  parameterization  has  some  advantages  in 
statistical  work  since  Gamma(k,y)  =  y  Gamma(k,l)  [Ref.  7]. 

A  Gamma (k,l)  random  variable  has  density 

V*’*'1'  -  rw  xk‘X  e'kX'  *  >  °  (II.B.1.2) 

This  is  a  special  case  of  the  more  general  density 

(k.k  kx 

fG(x;k,y)  =  TTkl"  x^_1  e~ 

where  y  =  1.  Of  course,  since  the  scale  parameter,  shape 

k 

parameter  and  mean  are  related  by  the  relation  y  =  the  den¬ 
sity  can  be  specified  by  any  two  of  the  three  parameters.  The 


typical  parameterization  in  terms  of  the  scale  parameter.  A, 


is  useful  because  GCk^A)  +  G(k2,A)  =  G(k1+k2,A).  This  re¬ 
lationship  is  not  true  when  the  parameterization  is  through 
the  shape  parameter,  k,  and  the  mean,  y. 

A  Beta(m,n)  random  variable  has  density 


fB (x;m,n) 


T (m+n)  m-1 

rTnOTTn)  x 


(1-x) 


n-1 


0  <  x  <  1,  m  >  0,  n>0. 


(II. B. 1.3) 


For  the  Beta  random  variable  to  be  properly  defined 

each  of  the  parameters  must  be  positive.  Hence,  when  q  =  k, 

(II. B. 1.1),  as  defined  above,  is  no  longer  appropriate  since 

each  Beta  random  variable  has  a  parameter  that  is  identically 

zero.  In  this  case  when  q  *  k  it  is  understood  that  the  {An> 

sequence  is  considered  to  be  identically  zero  and  the  {Bn} 

sequence  one.  Therefore,  I I. B. 1.1.  becomes  simply  X  =  G  , 

n  n 

and  the  (Xn)  sequence,  like  the  {Gn>  sequence,  is  iid.  A 
justification  of  this  generation  scheme  for  a  Gamma  process 
as  defined  by  II. B. 1.1  was  provided  by  Lawrance  and  Lewis 
[Ref.  14 ,  pp . 24 ] . 

In  this  section  the  correlation  structures  of  the  \ X  ’ 
sequence  and  that  of  the  {X^}  and  (G  ;  sequences  are  addressed. 
Other  characteristics  of  the  sequence,  such  as  conditional 
expectation  of  Xn  given  Xn_^,  directional  moments,  and 
P(Xn+1  '■  X  )  are  considered.  Of  particular  note  is  the  deri¬ 
vation  of  the  conditional  density  of  Xn  given  Xn_^.  This 


leads  to  the  formulation  of  the  likelihood  function  and  a 
computer  program  to  generate  maximum  likelihood  estimates 
of  parameter  values.  The  numerical  convergence  properties 
of  the  likelihood  method  are  assessed  in  Section  II. E. 

2 .  Correlation  Structure 

The  serial  correlation  of  the  { X^}  series  is  easily 
determined  by  a  straightforward  calculation.  We  have 


Now  Xi  and  G ^  are  independent  if  j  >  i  and  Xi  and  A^  are 
independent  if  j  >  i.  Using  these  facts  along  with  the  iid 
nature  and  independence  of  the  {Bn}  and  {Gn>  sequences  yields 
the  following  expression  when  expectations  are  taken. 


E  (A  )  E  (X^  .)  +  E  (B  )  E  (G  )  E  (X  ,) 
n  n-1  n  n  n-l 

<*j2>  •  +  (2)  -1-1 


Therefore , 


COV(Xr,,Xr,  ,) 
n  n-i 


and 
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CORR (X  , X  ) 
n  n-1 


0  <  q  <_  k. 


(II. B. 2.1) 


This  is  consistent  with  the  fact  that  for  q  =  k,  { X  }  is  a 

n 

sequence  of  i.i.d.  Gamma  variates  which  implies  that 

CORR (XX  .)  =0.  This  correlation  is  easly  extended  by  an 
n  n-i 

induction  argument  to  yield 


CORR (X  , X  ) 
n  n+m 


k  k  ’  ' 


n  >  m  >  0 . 


(II.B.2.2) 


The  two  sequences  {X  }  and  {G  }  can  be  viewed  as  a 

n  n 

bivariate  pair  (X  ,G  )  of  Gamma (k,l)  random  variables. 
Therefore,  the  correlation  structure  of  these  sequences  may 
be  of  interest.  Proceeding  as  before 


AX  .  +  BG 
n  n-1  n  n 


X  G 
n  n 


AX  , G  +  BnG„ 
n  n-l  n  n  n 


Taking  expectations  as  before 


E(X  G  ) 
n  n 


E  (A  )  E  (X  .  )  E  ( G  )  +  E  (B  )  E  (G^) 
n  n-1  n  n  n 


*<W  -  + 


E (xnG  ) 
n  n 


k  +q 

=  if 


(II.B.2.3) 


Therefore 


COV(X  ,G  )  *  -% 

n  n  k2 


and 
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(II.B.2.4) 


COER(Xn,Gn)  =  |  . 


When  q  =  k  this  is  1  since  X  =  G  . 

n  n 


Pursuing  the  process  one  more  step  we  determine 


CORK (X  , G  .). 
n  n-i 


X„  =  A  X  ,  +  B  G  , 
n  n  n-l  n  n 


X  G  -i  =  AX  -G  ,  +  B  G  G  , 
n  n-l  n  n-l  n-l  n  n  n-l 


Taking  expectations  as  before  and  using  II.B.2.3  and  the  second- 
order  stationarity  of  the  {X^}  sequence,  we  get 


E (X  G  .) 
n  n— l 


=  E(A  )EU  1  G  .)  +  E  (B  )  E  (G  )  E  (G  .) 
n  n-l  n-l  n  n  n-l 


Therefore, 


»  (IsiS)  (1l+3)  =  k_+ka-k.Vs2+qk2 

k  k2  k  k3 

k3+kq-q2 

k3 


COV(X  ,G„  ,) 
n  n-l 


q  (k-q) 


(II.B.2.5) 


and 


CORR (X  , G  .)  = 

n  n-l 


(II.B.2.6) 


II.B.2.3  through  II.B.2.6  can  be  used  in  a  simple  induction 
argument  to  yield  the  general  result 


CORR(X  , G  J 
n  n-m 


<|)  (^2)m'  m  =  O'1'  «  •  •  'n' 


(II.B.2.7) 


When  j  is  greater  than  i,  X^  and  G  j  are  independent.  Hence, 
CORR(Xi,Gj)  *  0,  j  >  i. 
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3 .  Conditional  Expectation  and  Conditional  Density 

The  conditional  expectation  of  X  given  X  ,  =  y  is 

n  n— x 

trivially  determined  from  trie  defining  equation  and  the 
moments  of  the  Beta  distribution  as 

E(XnlXn-l  =  y)  =  (^T2')Y  +  k*  (II. B. 3.1) 

Recognizing  ^3  as  the  correlation  between  Xn  and  Xn-1  and 
letting  p  be  that  correlation,  II. B. 3.1  can  be  written  as 

E(xn|xn_i  =  y)  =  py  +  |  =  py  +  (l-p).  (ii.g.3.2) 

Thus  the  regression  is  linear  in  y. 

The  conditional  density  of  Xn  given  Xn_^  can  also  be 
determined.  It  is  easiest  to  start  by  deriving  the  condi¬ 
tional  distribution  of  X  given  X  ,  =  y. 

n  n-i 

P  (Xn  <  x\x  ,  =  y)  =  P  (A  y  +  B  G  <  X)  =  P  (B  G  <  x  -  A  y) 
n  —  1  n-l  n  nn—  nn—  n 

Now  [Ref.  14]  the  product  of  a  Beta(q,k-q)  random  variable  and 
a  Gamma (k,l)  random  variable  is  a  random  variable  with  density 
Gamma(q,^).  Hence,  if  we  let  Dn  be  a  Gamma(q,^)  random 
variable, 


P (X  <  x  j  X  .  =  y)  =  P (D„  <  x-Av)  (II.B.3.3) 

n  —  '  n-l  n  —  n 

This  can  be  written  as  a  convolution  if  one  is  careful  about 

the  upper  limit  of  integration.  Since  D  is  Gamma(q,^),  it 

n 

is  non-negative.  Hence,  P(Dn£X-Any)  is  zero  if  x~Any  is 
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less  than  zero.  Since  An  is  a  Beta  random  variable  and, 
hence,  bounded  above  by  one,  x-Any  can  not  be  negative  if 
x  >  y.  However,  if  y  >  x,  then  Afl  must  be  restricted  to  lie 
in  the  range  (0,^) .  Taking  this  restriction  into  account 
and  writing  the  RHS  of  II.B.3.3  as  a  convolution 


p(Xn <  xj 


Xn-1  =  ^ 


L 

/  f  (z)F  (x-yz)dz, 
0 


and 


L 


1  if  x  ^  y, 
x  .  r 

y  if  x  <  y. 


(II.B.3.4) 


where  f A(z)  is  the  density  of  a  Beta  random  variable  as  in 
II. B. 1.3,  and  is  the  distribution  function  of  a  Gamma 
(q,^)  random  variable. 

Of  course,  to  get  the  conditional  density  of  XR 
given  Xn_^  =  y,  we  must  take  the  derivative  of  II.B.3.4 
with  respect  to  x.  Recognizing  that  the  upper  limit  may 
be  a  function  of  x  and  applying  Leibneitz's  Rule  where 
appropriate  yields 


(x) 

‘n-1 


L 


/ 

0 


fA(z)  f D (x-yz) dy , 


and 


L 


if  x  .>  y, 
if  x  <  y, 


(II .3.3.5) 
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i 


where  fv  ,v  (x)  is  the  conditional  density  of  X  given 
n '  n-1  n 

Xn_lf  the  Beta  density,  and  fD(x-yz)  is  the  Gamma 

(<3/^)  density  as  in  II. B.  1.2.  Writing  this  result  in  terms 

of  the  densities  involved  we  have 


fX  | X  .  {x) 
n 1  n-1 


=  fL  _ LljO  _ z^-q-l/1-  =  )q~1  kq  (X- v~)q~1n~k(x~yz) dv 

f  (x-q)  T  (q)  z  U  z)  T(a)U  y2)  e  dy ' 


with  the  condition  on  L,  as  in  II. B. 3.5.  And  finally, 


fX  | X  . (x) 
n 1  n-1 


F (k-q)  [ F  < q) ] 


5)k9e-kx 


.  ( 


x  /  zk"q~1(l-z)'3':l'(x-yz)q_1akyzdy,  (II.b.3.6) 

0 


1  if  x  >_  y, 


U  i 


f  x  <  y. 


As  a  check  on  the  derivation  of  the  conditional  density, 
the  conditional  density  and  conditional  expectation  were  calcu¬ 
lated  for  values  of  k  and  q  which  produced  simple  integrands. 


One  of  these  cases  was  2,  q  =  1.  Then  k-q-1  =  0  and 

q-1  =0.  In  these ' parameter  values  II.B.3.6  reduces  to 


£x  ,y  (x)  =  2e-2x  /L  e2yzdz  L  =  j 

Xn'Xn-l  0  ^ 


1  if  x  _>  y. 


^  —  if  x  < 

y 
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After  integration. 


-2x  re^  1,  .  , 

e  [  — - -]  if  x  >  y, 


fv  |  (x) 

xn '  n- 1  I  1  e-2x 

y  y~ 


(II.B.3.7) 


if  x  <  y. 


Since  the  two  expressions  in  II.B.3.7  are  non- negative,  we 
can  insure  that  this  is  a  density  by  verifying  that  its  inte¬ 
gral  is  one . 


/  fy  | „  (x) dx 


0  'XnlXn-l 


0  y  Y  y  y  Y 


=  -  dx  -  -  /7  e”2xdx  +  [ 

\7  J  V  ' 


.2y 


hi  e*2xdx 

y  y 


1  +  .  .L  *  _L  . 

1  2y  2y  2y  y 


=  1, 


as  required.  We  can  also  take  the  conditional  expectation  to 
see  if  it  equals  ^  +  j  as  required  by  II. B. 3.1.  Thus, 


/  xfv  |  v  (x)dx 

0  n '  n-l 


x[i- 

y 


e"2*,  .  ^  -2x.e2y 

— ]dx  +  /  xe 


-  — ]  dx 

y 


/  Xfv  lY 
0  xn|xn-l 


1  rY  lr*  -2x  e^y  1  ,  _2  v 

(x) dx  =  i  /  xdx  -  i/  xe  2xdx  +  [^—  -  ±]  /  xe  2xdx 


y  y 


y  ,  e'2y  .  e'2y 
2  +  2  +  4y 

e~2y  e-2y 


_L  +  I  +  JL 

4y  2  4y 


=  £  +  i. 
2  2' 


as  required. 


Using  k  =  3  and  q  =  1  produces  a  density  of 


fXlX  <=°  = 

n '  n-i 


_  -3x rl ,  3y  eJy  ^  1  ,  . e 

2e  [y(e  3y-  +  3$)]  lf  x  i  y  r 


2x  _2  n  ~~3x. 

-  7-7(1  -  e  ) 
y  3y 


(II.B.3.8) 


if  x  <  y  . 


This  density  is  non-negative  and  integrates  to  one.  It  also 
produces  a  conditional  expectation  of  ^  as  required  by 
II. B. 3.1.  These  results  are  also  of  use  in  validating  the 
results  obtained  from  numerical  integrations  of  II.B.3.6  in 
estimation  applications. 

This  conditional  density  can  be  used  to  form  the  joint 
density  of  . . .  ,XR  and,  hence,  the  likelihood  function 

of  the  data.  This  subject  will  be  addressed  in  the  following 
section. 

4 .  Maximum  Likelihood  Estimation 

Once  the  conditional  density  of  X  given  X  .  and  the 

n  n-l 

marginal  density  of  are  known,  it  is  possible  to  evaluate 
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f 


the  joint  density  of  ,x2 , • • • iXfi .  Since  the  first-order 
autoregressive  process  is  Markovian,  as  can  be  seen  directly 
from  the  defining  equation  II. B. 1.1,  the  following  equation 
is  valid: 


f  (  Xj  |  x j  _ ^ ,  x j  f  *  *  *  i  x^ ) 


f  (  Xj  I  Xj^)  ,  n  >_  j  >_  2  (II. B. 4.1) 


Applying  II. B. 4.1  n-1  times  to  the  joint  density  of  Xn,Xn_1» 
...,X^  produces 


f  (^^n-l'-**'^  =  fl(xnlVl)fl(Vllxn-2)  . . .  f  L  (x2 1  x)  f  2  (x^ 


(II. B. 4 .2) 


where  f  is  the  joint  density  of  Xn»...,X^;  is  the  condi¬ 
tional  density  of  X.  given  X  and  f  is  the  marginal  den- 

J  J  "" 

sity  of  the  (Xn)  sequence. 

Viewing  the  joint  density  as  the  likelihood  function, 
letting  L  be  the  likelihood  function,  and  taking  natural 
logarithms  of  each  side  of  II.B.4.2  produces 


n-1 

In  L  -  Inf,  (x,  )  +  l  In  f,(xu1|xJ  (II.B.4.3) 

*  i=l  -  1 

Recall  that  in  Section  II. B. 3  the  conditional  density 
was  determined  to  be 
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fjCxiy)  =  { 


r  (k) 


r(k-q)  tr(q) ] 


/  2^-'-(1-z)'3-l(x.yz)‘3-18ky2d2i 


,  ( 


1  if  x  21  y< 


£  if  x  <  y. 


(II. B. 4 .4) 


For  a  given  set  of  data  the  likelihood  can  be  viewed  as  a 
function  of  the  parameters  k  and  q.  Let 


G  (k,  q)  =  In  f-( 


n-1 


x  )  +  £  In  f.( 

1  i=l  1 


xi+l ' xi} 


(II. B. 4 .5) 


Assume  for  the  moment  that  a  procedure  has  been  established 
to  evaluate  G(k,a)  given  k,q  (0  <  q  £  k)  and  the  data.  Consider 
the  problem  of  constructing  a  program  to  determine  the  values 
of  k  and  q  which  maximize  G(k,q) .  An  outline  of  the  program 
can  be  constructed  as  follows: 

1.  Read  the  initial  values  for  k  and  q. 

Read  the  data. 

2.  Determine  a  direction  of  search.  Use  the  following 
difference  equations  to  approximate  derivatives. 


G ( k, q) 


G(k+Ak,q)  -  G(k,q) 
Ak 


(II. B. 4 .6) 


oq 


G  (k,  q) 


G  ( k , g+Aq)  -  G(k,q) 

Aq 


(II. B. 4 .7) 
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Let  Dk.  =  GjkjSL  and  Dq.  =  G(^>^.q)  -  G(k,g) 

1  Ak  Aq 

Dk 

for  the  ith  iteration.  If  we  define  7G^  =  (^i)  >  then  7G^ 
approximates  the  gradient  of  G  at  the  current  k,q  values. 

For  i  =  1,  let  the  direction  of  search,  d^,  be  7G^.  For 
i  >  1,  define  the  direction  of  search,  d^,  to  be 

dk  =  7Gi  +  Bk-1  dk-l'  (II.B.4.8) 

7GT7G. 

where  0,  1  is  — ^ — - ,  the  ratio  of  the  length  of  the 

7Gi-i7Gi-i 

present  and  preceding  gradients.  Formula  II.B.4.8  is  the  key 
equation  in  implementing  the  Fletcher-Reeves  Conjugate  Gradi¬ 
ent  Algorithm.  Once  d^  has  been  selected,  normalize  its 
length  to  one. 

_3 

3.  Let  the  initial  step  length,  SL,  be  10  and  let  N  =  1. 
Compute  the  trial  values  of  k  and  q,  Tk  and  Tq,  using 


Tk  =  k  +  SL  *  Dk. , 
Tq  =  a  +  SL  *  Dq^. 


(II.B.4.9) 


If  G (Tk, Tq)  >  G(k,q) ,  let  SL  =  2  *  SL ,  otherwise  go  to  4. 
If  N  =  10;  k  =  Tk,  q  =  Tq,  go  to  2.  (Ho  step  larger  than 
210  *  10-3=1.0  is  allowed.)  If  n  <  10;  N  =  N+l,  go  to  3. 

4.  If  N  >  1  (  at  least  one  step  produced  an  increase), 
use  a  golden  section  search  between  Tk,Tq  for  step  N-2  and 
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Tk  and  Tq  for  step  N  to  determine  the  maximum  function 
value  and  the  k  and  q  values,  kMAX  and  qMAX,  which  produced 
it.  Here  k  =  kMAX,  q  =  qflAX,  go  to  2. 

If  N  =  1,  go  to  5 . 

5.  Since  the  initial  step  along  the  indicated  direction 

produced  a  decrease  in  the  function  value,  check  to  see  if 

you  are  at  a  local  maximum.  Determine  the  function  value 

at  points  at  30°  intervals  on  the  circumference  of  circles 

-3  -2  -1 

with  radii  of  10  ,10  ,  10  (0°  is  parallel  to  the  q  axis) . 

If  the  maximum  of  these  test  values  is  greater  than  the  pre¬ 
sent  value,  set  k  and  q  to  the  values  which  produced  the 
maximum  value,  set  the  present  function  value  to  the  maxi¬ 
mum,  set  i  =  1  and  go  to  2.  Otherwise  terminate. 

The  program  above  assumed  that,  given  q,  k  and  the 
data,  the  value  of  the  likelihood  could  be  calculated.  One 
difficulty  in  performing  the  calculation  is  that  the  integrand 
of  the  conditional  density  may  contain  singularities.  As  a 
precondition  to  using  an  IMSL  routine  to  evaluate  the  inte¬ 
gral,  these  potential  singularities  must  be  removed.  The 
technique  employed  requires  that  the  coefficient  of  the  term 
that  goes  to  infinity  as  one  of  the  limits  of  integration  is 
approached  is  added  and  subtracted.  The  part  that  is  added 
is  then  integrated  separately  and  added  to  the  part  that  is 
evaluated  by  the  IMSL  routine. 

To  procede  with  this  technique  we  first  split  the 
integral  into  two  parts.  Thus,  ignoring  the  part  of  II.B.4.4 
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outside  the  integral,  we  have 


yz^'V^dz  - 

0  0 

x  (x-yz)q  ^ekyzdz  +  /  zk  q  ^(l-z)q  ^(x-yz)q  ^ekyzdz 

L/2 

(II. B. 4. 10) 

In  the  first  part  of  the  RHS  of  II. B. 4. 10  the  term  z  H  could 

tend  to  infinity  as  z  tends  to  zero  if  k-q-1  <0.  If  we  set 

z  equal  to  zero  in  the  remainder  of  the  integrand  we  get 

(1-z) q_1 (x-yz) q-1ekyz |  =  xq_1.  Adding  and  subtacting  this 

z=0 

term  times  the  term  that  contains  the  singularity,  we  have 


L/2 

/  zk-'!-1(l-zlq-1(x-yz)'’-1elt''zdz 

0 


=  jL/2zk-q-1(l-z)q-1(x-yz)q-1ekyz-!<q-1z,c-q-Lx<3-12k-q-1dz 

0 


=  (L/2zk~q-i  [  (l-z)q-i  (x-yz)q-1ekyz-xq-1']  dz  +  xq’1  {L/2zk'q-ld 
0  0 

Thus , 

L/2 

/  zk_q_1(l-z)q'1(x-yz) q_1ekyzdz 

0 


.  /L/2zk-q-1[(l-2)q'1(x-yz)q-1ek!'z-xq-1ldz«q-1 

0 


(l)k-q 

k-q 


(II. B. 4. 11) 


Recall  that  since  q  <  k,  k-q  >  0  and  the  second  integral  on 
the  RHS  can,  in  fact,  be  integrated  as  shown.  Now  the  inte¬ 
gral  on  the  RHS  can  be  evaluated  by  IMSL  routine  DCADRE  and 
the  second  part  of  the  RHS  can  be  easily  computed. 

Applying  this  technique  to  the  second  half  of  the 
integral,  recalling  that  the  case  where  L  =  1  and  L  =  — 
must  be  considered  separately  produces 

f1  zk-<3-1(l-z)'3-1(X-y2)<J-1e!t5'2dZ 
1/2 


==  /"  (l-zl^az  (II.B.4 

1/2 

g-i  ky  <l/2>q 
+  (x-y)q  eKy  q 


and 

jX/Y  z1'-’-1  (1-z)  I-1  (x-yz)<J-1ekyzdz 
x/2y 


=  /X/Y  (x-yz)<3-1[zil-<5-1(l-z)c'-1ek',z-(?)k'‘3':1(l-S)<J'1ekx]dz 
x/2y  y  y 


+  (— )  k-q_l  ( i_X)  q— 3-ekx 


(|)q 

yq 


(II.B.4 


.12) 


.13) 


Since  q  >  0,  the  second  terms  in  the  two  previous  expressions 
are  properly  integrated. 


Two  points  should  be  noted.  First,  if  k-q-1  >  0  or 
if  q-1  >  0,  then  these  steps  are  not  required.  But  whether 
they  are  required  or  not,  they  are  always  accurate.  Since 
the  exact  path  of  the  search  algorithm  is  unknown  at  the 
start,  these  expressions  are  used  throughout  to  insure  accu¬ 
rate  calculations  regardless  of  the  k  and  q  values  encoun¬ 
tered.  Second,  if  x  =  y,  two  parts  of  the  integrand  simul¬ 
taneously  tend  to  infinity  and  this  procedure  breaks  down. 

This  does  not  pose  a  problem  for  continuous  data  since  the 
probability  of  this  occurring  is  zero.  However,  if  discrete 
values  are  used  or  if  the  data  is  truncated  because  of  limits 
on  the  accuracy  of  the  measuring  device,  then  the  data  may 
have  to  be  preprocessed  to  insure  that  successive  values  are 
not  equal  by  adding  a  small  increment  to  one  of  the  values. 

When  the  program  was  written,  its  accuracy  was  veri¬ 
fied  by  three  checkes .  The  case  k  =  q  implies  independence 
in  the  basic  model  since  as  q  tends  to  k  the  probability  that 
An  equals  zero  tends  to  one  and  the  probability  that  Bn 

equals  one  tends  to  one.  Hence,  in  the  limit  X  =  G  and  G 

n  n  n 

is  an  iid  sequence.  The  logarithm  of  the  likelihood  function 
for  independence  was  calculated  and  compared  to  the  program 
results  for  several  values  of  k  and  q  where  k  =  q.  The  two 
calculations  were  equal  within  machine  roundoff  and  compu¬ 
tational  accuracy.  The  special  cases  of  k  =  2,  q  =  1  and 
k  =  3,  q  =  1  discussed  in  II. B. 3  were  also  computed.  The 
logarithm  of  the  likelihood  function  was  computed  using  each 


of  the  conditional  densities  derived  in  II.B.3.7  and  II.B.3.8. 
When  the  results  of  these  calculations  were  compared  to  the 
program  results  with  the  specified  k  and  q  values,  the  re¬ 
sults  were  equal  within  calculation  and  roundoff  accuracy. 

The  results  of  the  tests  of  this  program  when  used 
with  simulated  data  with  known  parameter  values  are  presented 
in  Section  II. E. 

Note  that  there  are  natural  moment  estimators  for  the 
three  parameters,  k  and  q  and  •„  in  this  model.  These  follow 
from  the  fact  that 


0(1)  =  CORR(XnXn+1)  •  1  -  f, 

c2(x)  =  Var(Xn)  =  var(x)  =  °X  =  1 

[E(Xn)]2  [E (X) ] 2  k 

Thus  we  use  for  moment  estimations 


u  =  x 

q  =  (l-p(l)k 


(II. B. 4. 14) 

(II. B. 4. 15) 

(II .3.4.16) 


These  moment  estimations  will  be  compared  to  maximum  likeli¬ 
hood  estimations  in  the  case  where  y  =  1  in  Section  II. E. 


5 .  Directional  Moments 

Unlike  processes  with  normal  marginals,  non-normal 
processes  are  not  completely  determined  by  their  correlation 
structure.  Directional  moments  not  only  demonstrate  this 
difference,  but  also  help  to  differentiate  processes  with 
similar  correlation  structure  and  identical  marginal  distri¬ 
butions.  They  can  also  be  used  to  help  determine  parameter 
values.  In  addition,  they  may  also  be  viewed  as  another 
way  of  characterizing  the  joint  distribution  of  the  process. 
With 


=  A  X 


n  n-1 


+  BnG. 


n' 


with  all  random  variables  defined  as  in  II. B. 1.1,  we  first 
2 

address  E(X  X  . )  .  We  have 
n  n— l 


Xn  =  A  X  .  +  B  G 
n  n  n— 1  n  n 


XX2.  =  A  X3  +  B  G  X2  , 
n  n-1  n  n-1  n  n  n-1 


Taking  expectations,  recalling  G^  and  X^  are  independent  if 
i  >  j,  { Bn }  and  { Gn }  are  independent,  and  A^  and  Xj  are 
independent  if  i  >  j .  Then 


E(Vn-l> 


E  { A  )  E (X3  .)  +  E (B  ) E (G  ) E (X2  .) 
n  n— 1  n  n  n-1 


k-q  k (k+1) (k+2)  +  ^ 


k 


k (k+1) 
,2 


1 


(k+2)  +  q] 


E(Wn-l> 


_  k (k+1)  r (k-g) 

~  3  1  k 

kJ  * 


E(X  X„ 
n  n-1 


k (k+1)  rk  +2 (k-q) , 

-  - 5 — l - r — 

kJ  K 


(II. B. 5.1) 


Since 


2  2  2  2  2 
X„  =  .  +  2ABGX  ,  +  B„G* 

n  n  n-l  n  n  n  n-I  n  n 


2  2  3  2  2  2 

X„X„  .  =  A„X„  .  +  2A  B  G  X^  ,  +  B„G^X„  , 

n  n-l  n  n-l  n  n  n  n-l  n  n  n-l 


Taking  expectations  as  before 


E (X^X  .)  *  E(A^)E(X^  .)  +  2E ( A  ) E  (B  ) E (G  ) E (X^  .) 

n  n-l  n  n-l  n  n  n  n— I 

+  E(B^)E(G^)E(Xn_1) 


(k-q) (k-q+1)  k (k+1) (k+2)  ,  2 (k-q)  q  ,  k(k+l) 

k (k+1)  k3  k  'k'1'  k2 


+  q(q+1.)  Mkjll.i 

k(k-l)  k2 

After  simplification  we  get 

E(XnXn-l)  =  (k+2)  +  (II.B.5.2) 

k3  k3 

Note  that  these  two  directional  moments  are  different, 
indicating  that  the  process  is  not  time  reversible. 
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6-  P(xn+l>xn> 


Another  characterization  of  the  joint  distribution 
and  the  time  directionality  of  the  process  is  given  by 
P (Xn+^  >  X  ) .  There  is  a  simple  analytical  solution  for 
P  (Xn+^  >  Xn)  in  the  GLAR(l)  process.  Consider, 


P!xn+l>Xn>  " 


p(Vlxn 


B  ,,  G  >  X  ) 
n+1  n+1  n 


=  p(VlCn+l>  [1-An+l1X„»  (II.B.6.1) 

Recall  that  Bn+^  is  Beta(q,k-q)  and  An+^  is  Beta(k-q,q)  . 

Hence,  [1-An+^]  is  Beta(q,k-q)  .  Since  Gn+^  and  Xn  are  inde¬ 
pendent  and  have  the  same  marginal  distribution  and  An+1  and 
Bn+1  are  independent,  each  side  of  the  inequality  in  XI. B. 6.1 
has  the  same  distribution  and  the  random  variables  are  inde¬ 
pendent.  Hence 


p<Vi>x„>  -  °-5°- 

This  is  a  strong  property  of  the  process.  While  the 

process,  as  seen  by  the  directional  moments,  is  not  time 

reversible,  the  fact  that  P(X  ,,  >  X  )  =  0.50  means  that 

n+i  n 

sample  paths  will  have  as  many  "runs  up"  as  "runs  down". 
Sample  paths  are  given  in  Figures  II.B.6.1  and  II.B.6.2. 

An  additional  point  of  interest  occurs  when  k  =  1 
and  the  process  has  exponential  marginal  distributions. 


4  7 
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Another  exponential  first-order  autoregressive  process  in 
which  p(xn+i>xn)  =  0.50  is  the  PREAR(l)  process.  This  is 
a  special  case  of  the  two  parameter  NEAR ( 1 )  exponential 
process  of  Lawrance  and  Lewis  [Ref.  8  ]  in  which  the  two 
parameters  a  and  8  are  related  by  8  =  .  The  two 

exponential  processes  are  very  similar  in  sample  path  proper¬ 
ties.  However,  the  GLAR(l)  process  has  a  smoother  joint 
distribution.  In  fact,  the  likelihood  for  the  PREAR(l) 
process  is  discontinuous,  making  it  difficult  to  get  maximum 
likelihood  estimators. 

C.  FIRST-ORDER  MOVING  AVERAGE  BETA-GAMMA  MODEL,  GLMA(l) 

1 .  Introduction 

Another  special  case  of  the  GLARMACp , q)  model  is  the 
first-order  moving  average  model  where  p  *  0  and  q  =  1.  This 
arises  naturally  from  the  key  result  that  an  { Xn }  sequence 
can  be  formed  by  the  random  sum  of  two  independent  Gamma  ran¬ 
dom  variables.  In  the  first-order  autoregressive  case  of 
Section  II. B,  the  generation  scheme  was  given  by  equation 
II. B. 1.1  and  is  repeated  here 


n 


=  AX  .  +  B  G  . 
n  n-1  n  n 


The  distribution  of  the  {x  }  sequence  depends  on  the  inde¬ 
pendence  and  distribution  characteristics  of  X  ,  and  G  . 

n-x  n 

It  should  be  noted,  however,  that  any  two  independent  random 
variables  with  the  required  Gamma  distribution  could  be 
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- 


substituted  for  X  ,  and  G  without  changing  the  distribution 

n-i  n 

of  X  .  In  particular  if  we  substitute  G  for  X  ,  and  G  , 
n  c  n  n-1  n-1 

for  Gn,  we  produce  the  first-order  moving  average  process 
which  generates  the  {X^}  sequence  using 


xn  *  AnGn  +  BnGn-l'  (II. C. 1.1) 

where  {Xn,  n  =  1,2,...}  is  a  second-order  stationary  sequence 

of  random  variables  with  marginal  Gamma{k,l)  distributions, 

{A^,  n  =  1,2,...}  is  an  iid  sequence  of  Beta (k-q,q)  random 

variables;  (B  ,  n  =  1,2,...}  is  an  iid  sequence  of  Beta (q,k-q) 

random  variables;  (Gn,  n  =  0,1,...}  is  an  iid  sequence  of 

Gamma  (k,l)  random  variables;  {A  },  (B  } ,  and  (G  }  are  inde- 

n  n  n 

pendent;  0  <  q  £  k.  The  Gamma  random  variables  are  para¬ 
meterized  as  in  II.B.l  with  density  as  in  II. B. 1.2.  The 
Beta  random  variables  have  density  as  in  II.B.l. 3. 

In  this  section  we  will  address  the  correlation  struc¬ 
ture  of  the  {X  }  sequences  and  that  of  the  (X  },  (G  }  se- 
n  ^  n  n 

quences,  theoretical  ranges  for  possible  correlations  for 

the  {X  }  sequences,  directional  moments,  and  the  P (X  +1  >  Xn) . 

2 .  Correlation  Structure 

Using  II. C. 1.1  to  define  X  and  X  ,  we  have 
3  n  n-i 


XX  .  =  (AG  +  B  G  ,)  (A_  ,G_  ,  +  B„  .  G_  7) 

n  n— i  n  n  n  n-i  n-i  n— i  n— i  n- z 


=  AA  .  G  G_  i  +  A  B  .  G  G  0  +  A  .  B  G  ,+BB  ,G  ,G  0 , 
n  n-1  n  n-l  n  n-1  n  n-2  n-1  n  n-1  n  n-1  n-1  n-2 
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Using  the  iid  nature  and  independence  of  {A  },  {B  },  and 

n  n 

{ Gn }  and  taking  expectations 

E<xnxn-1>  -  E(VE(An-l,E(Gn)E(Gn-l)+E(VE<Vl>E(Gn)E(Gn-2> 

+  ElAn-l)E(Bn,E(Gn-l,+E(Bn)E,Bn-l)E(G„-l)E(Gn-2) 

-  't3'  2  +  't2'  <?>  +  (t2'  <f>  +  (?) 2 

k 

-  1  +  (^)  (|)  <£>  (II.  C.  2.1) 

Therefore, 

C0V(VVl>  =  It3)'?)1!) 

and 


CORR(Xn,Xn_1)  =  <1 -§)(§>•  (II.C.2.2) 

A  simple  calculation  will  show  that  for  lags  greater 
than  one  the  correlation  is  zero.  So  equation  II.C.2.2  plus 
the  knowledge  that  greater  lags  are  zero  is  sufficient  to 
specify  the  correlation  structure  of  the  {Xn}  sequence. 

One  might  note  at  this  juncture  that  the  correlation 
of  the  (X  }  sequence  is  constrained  to  lie  in  the  interval 
(0,j)  .  The  reason  for  this  constraint  and  a  method  for  re¬ 
laxing  it  will  be  discussed  later.  It  is  also  worthy  of 
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note  that  the  {Xn>  sequence  is  stationary  and  has  the  same 

marginal  distribution  as  the  {Gnl  sequence,  if  XQ  =  GQ. 

As  indicated  in  II. B. 2  the  {X  }  and  {G  }  sequences 

n  n 

can  be  considered  to  be  a  bivariate,  correlated  Gamma(k,l) 
process.  As  such,  the  correlation  structure  of  this  bivari¬ 
ate  Gamma  may  be  of  interest.  Consequently,  we  first  develop 

the  correlation  of  X  and  G  in  the  standard  fashion. 

n  n 


n 


AG  +  B  G  t 
n  n  n  n-1 


X  G 
n  n 


A  G  +  B  G  G  , 
n  n  n  n  n-1 


Taking  expectations  as  before, 


E(W  -  E(VE<Gn>  +  E(BnIE(GniE(Gn-l) 

=  <^2>  ciUKJiX)  +  | 

k^  +  k  -  o 
k2 

Therefore, 


COV (X  , G  ) 
n  n 


and 
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(II.C.2.3) 


CORR(Xn,Gn) 


Now  consider  the  correlation  of  X  and  G  . 

n  n— l 


We  have 


X„  =  AG  +  B  G  , 
n  n  n  n  n-1 


X  G  . 
n  n— i 


A  G  G  ,  +  B  G  ,  , 
n  n  n-1  n  n-1 


Expectations  in  the  standard  fashion  produce 


E(XnVl>  =  E(An)E(Gn)E(Gn-l}  +  E(Bn)E(Gn-l> 


k  k  k2 


3if±SL 

k2 


Hence, 


COV (X  , G  ,)  = 

n  n— l 


.2- 

k2 


And  finally 


c°RE(Xn,Gn.i)  -  2 


(II.C.2.4) 


A  simple  calculation  convinces  one  that  the  correlation  for 

lags  greater  than  one  is  zero.  In  addition,  it  is  clear 

that  CORR (X  , G„,  )  =  0  for  m  =  1,2,...  Hence,  II.C.2.3  and 
n  n+m 


II.C.2.4  are  sufficient  to  specify  the  correlation  structure 
of  the  (Xn)  and  (Gn>  sequences. 

It  has  been  noted  before  that  the  range  of  correla¬ 
tions  for  the  first-order  moving  average  process  generated 
by  the  Beta-Gamma  method  of  II. C. 1.1  is  constrained  to  lie 
in  the  interval  from  zero  to  one-quarter.  There  may  be  other 
random  linear  combinations  of  Gamma  random  variables  which 
give  a  moving  average  process  with  Gamma  marginals  and  a 
greater  range  of  positive  correlations.  Thus  we  now  examine 
a  more  general  hypothetical  generation  process  to  prove  that 
any  random,  linear  combination  of  two  independent  Gamma  random 
variables  which  generates  a  sequence  with  the  same  first  two 
moments  as  those  Gamma  variables  has  a  correlation  that  lies 
in  this  same  interval.  In  fact,  this  proof  only  requires 
that  the  dependent  random  variable  have  the  same  first  two 
marginal  moments  as  the  generating  Gamma  variables. 

THEOREM: 

If  the  { Xn }  sequence  is  generated  by 

Xn  -  Vn  *  VW  (II.C.2.5) 

where  {Xn}  is  a  second-order  stationary,  non-negative  sequence 
of  random  variables  with  the  same  first  two  moments  as  the 
{G^}  sequence;  {An}  and  { }  are  iid  sequences  of  random 
variables;  {G  }  is  an  iid  sequence  of  Gamma (k,l!  random 


variables;  and  {A  } ,  {B  },  and  {G  }  are  independent,  then 

n  n  n 

0  <  CORR(X_ , X_  ,)  <  0.25. 

—  n  n-x  — 

Proof:  Since  {A  } ,  {B  } ,  and  {G  }  are  independent  and  {X  } 

-  n  n  n  n 

is  non-negative,  {An>  and  {Bn}  must  be  non-negative.  Hence 
E  (A)  _>  0  and  E(B)  >_  0 . 

Taking  expectations  of  II.C.2.5,  we  have 

E(Xn)  =  E(An)E(Gn)  +  E(Bn)E(Gn_1) 

1  =  E (A)  +  E (B)  (II.C.2.6) 


Hence,  0  <_  E(A)  <_  l  and  0  <_  E(B)  £  1. 

Computing  the  serial  correlation  of  1 Xn>  yields 


X  X  .  =  (AG  +  B  G  .)  (A  .C-  .  +  B  ,  Gn  ,) 

n  n— l  n  n  n  n-i  n— 1  n-l  n-i  n— 2 


—  A  A  .  G  G  -j  +  A  B  i  G  G  o  +  A  3  G  -i  +  B  B  .  G  G  ^ 
n  n-l  n  n-l  n  n-l  n  n-2  n— 1  n  n— 1  n  n-l  n-l  n-2 


Then 


£(X  X 
n  n- 


E  (An> E  lAn-l>  ElGn>  E(Gn-l>  +  E  (An> E  (Bn  -!»  E  <Gn> E  (Gn-2> 
+  E(A„-l)E(Bn)EIGn-l>  +E(Bn,E(Bn-l)E(Gn-l)E(Gn-2> 


=  1  +  E  (A)  E  (B) 
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Since  E(X  )  =  E(G  )  =  1, 

COV(Xn,Xn_1)  =  E(A)E(B)(|) 

And  since  VAR(X  )  =  VAR (G  ) , 
n  n 

CORK (X  , X  .)  =  E (A) E (B) 

n  n— l 

Using  II.C.2.6  and  its  consequences,  we  have 

0  <  CORR(Xn,Xn_1)  =  E (A)  [1  -  E (A)  ]  £  j.  (II.C.2.7) 

So,  in  general,  if  {X^}  is  second-order  stationary  with  the 
same  first  two  moments  as  {Gn>,  the  serial  correlation  of 
{Xn>  is  bounded  below  by  zero  and  above  by  one-fourth.  Q.E.D. 

This  constraint  on  the  correlation  appears  to  be 
restrictive  since,  in  the  classical  case,  when  two  normally 
distributed  random  variables  are  added  to  produce  a  normal  se¬ 
quence,  the  range  of  correlations  is  The  two  situa¬ 

tions,  however,  are  not  comparable.  It  is  clear  upon  reflection 
that  the  constraints  imposed  on  the  { Xn }  sequence  in  the  pre¬ 
vious  theorem  are  more  severe  than  those  imposed  upon  the 
classical  normal  case.  In  the  above  theorem  we  required  that 
both  the  mean  and  variance  of  the  {Xn }  sequence  equal  that  of 
the  innovative  sequence.  However,  in  the  classical  normal 
case  (where  zero  mean  normals  are  used  as  innovative  factors) 
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only  the  mean  of  the  generated  sequence  is  equal  to  the  mean 
of  the  innovative  sequence.  The  variances  are  not  equal. 

We  now  examine  the  case  where  the  generated  sequence  is  re¬ 
quired  to  have  the  same  mean  as  the  innovative  sequence,  but 
is  free  to  have  a  different  variance.  (This  is  the  case 
with  the  usual  constant  coefficient,  linear  additive  MA(1) 
scheme . ) 

THEOREM 

If  the  non-negative  (Xn>  sequence  is  generated  by 
II.C.2.5  and  all  variables  are  defined  as  for  that  equation 
except  that  { Xn }  is  only  constrained  to  have  the  same  first 
moment  as  {Gn},  then  0  £  CORR(Xn, Xn_^)  £  0.5. 

Proof :  Taking  expectations  of  II.C.2.5  with  the  new  circum¬ 
stances  produces 

E  (X)  =  [E  (A)  +  E(B)  ]  E  (G) 


1  =  E  (A)  +  E  (B) 

and  0  £  E(A)  £  1,  0  £  E (B)  £  1  by  following  reasoning  identi¬ 
cal  to  that  above.  Calculation  to  determine  the  serial  corre¬ 
lation  can  initially  proceed  as  usual. 

Vn-l  -  Vn*V«-lllVlVltVA-2l 


E(X  X  i> 
n  n-i 


1  +  E  (A)  E  (B)  (£) 
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COV(X  ,X„  ,) 
n  n-1 


=  E  (A)  E  (B)  (£)• 


To  this  point  all  calculations  and  reasoning  have  been  the 

same  as  that  which  produced  II.C.2.7.  However,  since  { Xn } 

is  not  constrained  to  have  the  same  second  moment  as  {G  } 

n 

the  most  explicit  result  obtainable  is 

E  (A)  £  (B)  (4) 

CORR(X„,Xn_i)  =  - VAR(X)—  '  (II.C.2.8) 


where  VAR(X)  is,  of  course,  a  function  of  VAR (A) ,  VARtB)  , 
and  VAR (G) .  Since  it  is  obvious  that  the  smaller  the  value 
for  VAR ( X ) ,  the  greater  the  serial  correlation  for  {Xnl, 
let  us  reduce  VAR(X)  to  its  smallest  values.  Since  the  dis¬ 
tribution  of  { Gn >  has  been  specified,  its  variance  is  fixed. 

Let  P (A  =  a)  =1  and  P(B  = b)  =1.  Then  trivially  E (A)  =  a, 
n  n 

E (B)  =  b  and  VAR(X)  =  (a^ + b^) (i) .  Under  these  conditions 
II.C.2.8  becomes 


CORR ( X, X  ,) 
n  n— i 


If  we  further  specify  that  a  =  b,  then  CORR(Xn,Xn_^)  achieves 

its  maximum  of  1/2.  Q.E.D. 

The  situation  developed  above  is  comparable  to  the 

classical  situation  where  the  innovative  factors  have  distri- 
2 

bution,  N(0,a  ).  Except  for  the  degenerate  case  where  one 
coefficient  is  zero,  the  sequence  generated  by  a  linear  com¬ 
bination  of  innovative  factors  may  have  the  same  first  moment 
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as  the  innovative  factors,  but  will  have  a  different  second 

moment.  So,  under  comparable  conditions  the  random,  linear 

combination  of  Gamma(k,l)  random  variables  can  produce 

positive  correlations  equal  in  magnitude  to  the  positive 

correlation  produced  by  a  classical  normal  process.  The 

distribution  of  the  {X  }  under  these  circumstances  is  unknown. 

n 

If  the  distribution  of  the  {X  }  is  constrained  to 

n 

be  Gamma (k,l)  and  the  Beta-Gamma  generation  scheme  is  used 
to  generate  the  {Xnl,  then  the  maximum  correlation  that  can 
be  achieved  is  one- fourth. 

3 .  Directional  Moments 

As  mentioned  in  II. B. 5  the  directional  moments  of  a 

non-normal  process  are  not  necessarily  equal  and  may  provide 

valuable  information  about  a  time  series.  First,  consider 

E(X2X„  .) .  From  II. C. 1.1 
n  n-i 


2  2  2  2  2 
X"  =  A„G„  +  2ABGG  .  +  B„G„  . 
n  n  n  n  n  n  n-l  n  n-l 


and 


X 


n-l 


A  .G  .  +  B  ■>  G  n 
n-l  n-l  n-l  n-2 


Therefore, 


X2Xn  .  =  A2An  ,G2Gn  -I  +2A  A  .  B  GG2  , +A  ,B2G2  , +A  B  , G*G 

n  n-l  n  n— l  n  n— l  n  n-l  n  n  n— l  n— l  n  n— l  n  n  i  n  n- 


+  2A  B  B  ,G  G  ,G  -+B2B  . Q2  . G  - 

n  n  n-l  n  n-l  n-2  n  n-l  n-l  n-2 


60 


Taking  expectations  as  in  II. C.  2  yields 


e(x2x  )  =  (k-q)  2  (k-q+1)  +  2(k-q)  2  (k+l|  +  (k-q)q(q+l)  (k+2) 


n  n-l 


+  (k-q)  (k-q+1)  q  +  2(k-g)2q2  +  q2(q+l) 
k3  k3  k3 


Upon  simplification  this  produces 


E(XnVl> 


.  ■^-[3k-q.3].1X~1>13[2'i(lI+1i.'|-k4'2l 


+  lr[2k-q+l] 


(II. C. 3.1) 


In  an  analogous  fashion 


xx2  .  =  A  A2  . G  G2  ,  +  2AA  , B  .  GG„  , G„  ,  +  A  B2  . G  G2  „ 

n  n— i  n  n-l  n  n— 1  n  n-l  n-l  n  n— 1  n-2  n  n— 1  n  n— 2 


+  A2  ,B  G3  ,  +  2A  B  B  .  G2  .  G  +  B  B2  .  G  ,G2  , 
n-l  n  n-l  n-l  n  n-l  n-l  n-2  n  n-l  n-l  n-2 


Taking  expectations  we  have 


E (X  X  .  ) 
n  n-l 


(k-q) 2 (k-q+1)  +  2 (k-q) 2q  +  (k-q) q (q-H) 
k3  k3  k3 


+  (k-q)  (k-q+2)  q  (k+2)  +  2  (k-q)  q2  (k+1)  +  q2  (q+1) 
k4  k4  k3 


which  simplifies  to 
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(k-q)  (k+q+1)  +  q(q+l)  (2k-q) 

k3  k3 


(k-q) q (k+1) 
k4 


(II.C.3.2 


4.  Empirical  P(Xn+1 >  X n> 

No  analytical  procedure  was  found  to  determine 
P(Xn+^>Xn).  Hence,  a  simple  computer  program  was  constructed 
to  evaluate  this  condition  for  a  series  of  sixty-eight  thou¬ 
sand  pairs  of  numbers  generated  by  the  Beta-Gamma  scheme  for 
each  of  ten  random  number  seeds.  The  answer  obtained  was 
considered  to  be  accurate  within  0.001.  The  comparisons  were 
run  for  each  of  seventy-nine  values  of  q,  from  0.05  to  3.25 
in  increments  of  0.05.  All  of  the  results  of  the  comparisons 
fell  in  the  range  0.499  to  0.501.  Fourteen  of  the  values  were 
different  from  0.500.  No  pattern  was  apparent  in  the  devia¬ 
tions  from  0.500  and  these  deviations  were  considered  to  be 
random  fluctuations  within  the  given  margin  for  error.  It 
thus  seems  clear  that  P(Xn+^>Xn)  for  this  process  is  like 
the  GLAR(l)  process  but  no  proof  has  been  found. 

D.  OTHER  CASES  OF  THE  GLARMA ( p , q )  MODEL 
1 .  Introduction 

A  primary  advantage  of  the  GLARMA (p,q)  model  is  the 
ease  with  which  it  can  be  adopted  to  cover  a  variety  of 
special  cases.  Two  special  cases,  the  first-order  autoregres¬ 
sive  GLAR(l)  and  the  first-order  moving  average  GLMA(l),  were 
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covered  in  II. B  and  II. C.  The  intention  here  is  to  briefly 
present  three  additional  cases  of  the  general  model  and  derive 
the  correlation  structure  of  each  case.  The  special  cases 
considered  are  the  first-order  mixed  model,  GARMA (1,1) , 
the  second  order  autoregressive  model  GAR(2) ,  and  a  bivariate, 
first-order,  autoregressive  model  BGAR(l) .  The  purpose 
in  presenting  these  cases  is  to  demonstrate  the  flexibility 
of  GLARMA ( p , q )  and  not  to  present  a  complete,  detailed  dis¬ 
cussion  of  each  model .  Further  extensions  of  the  special 
cases  of  the  GLARMA (p,q)  model  from  the  examples  given  are 
obvious.  Details  will  not  be  given. 

2.  GLARMA (1,1) 

Consider  the  following  scheme  for  generating  an  {Xn} 
sequence  of  random  variables. 


n 


—  BA  -,+CG, 
n  n-1  n  n' 


(II. D. 2.1) 


n 


=  DA  .  +  F  G  . 
n  n-l  n  n 


(II.D.2.2) 


where  {Xn,  n  =  1,2,...}  is  a  second-order  stationary  sequence 

of  Gamma (k,l)  random  variables;  {An,  n  =  0,1,...}  is  a 

second-order  stationary  sequence  of  Gamma  (k,l)  random 

variables;  (Bn,  n  =  1,2,...}  is  an  iid  sequence  of  Beta(k-q,q) 

random  variables;  {Cn,  n  =  1,2,...}  is  an  iid  sequence  of 

Beta  (q,k-q)  random  variables;  {D  ,  n  =  1,2,...}  is  an  iid  se- 

n 

quence  of  Beta(k-r,r)  random  variables;  {Fn, 


n  =  1,2,...}  is 


an  iid  sequence  of  Beta(r,k-r)  random  variables;  {B  },  {C  }, 

n  n 

{Dn>,  { Fn } ,  and  { }  are  mutually  independent;  0  <  q  _<  k; 

0  <  r  <_  k.  The  Beta  (m,n)  density  is  given  in  II. B.  1.3;  the 

Gamma(k,l)  density  in  II. B. 1.2. 

Before  the  serial  correlation  of  the  {Xn}  sequence 

can  be  determined,  the  serial  correlation  structure  of  the  {A  } 

n 

sequence  and  the  correlation  structure  between  the  {An>  and 
(Gn)  sequences  must  be  derived.  Proceding  first  with  the 
serial  correlation  of  the  {An>  sequence,  from  II.d.2.2  we 
have 


n 


=  DA  .  +  FG 
n  n-l  n  n 


So 


A  A  . 
n  n-l 


DA  .  +  FGA  . 
n  n-l  n  n  n— i 


Using  the  iid  nature  of  {D  },  {F  },  {G  };  noticing  that  when 

n  n  n 

i  >  j ,  D^  and  Aj ,  F^  and  Aj ,  and  G^  and  Aj  are  independent; 
recalling  the  independence  of  {Fn} ,  {Gn>;  and  taking  expec¬ 
tations  yields 


ElAnVl> 


E <A„A  1  > 

n  n-l 


,k-rx k (k+1)  r 

'  k  '  ,2  k  ' 

JC 

k2+k-r 


(II.D.2.3) 
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Therefore, 

C0V(VVl>  = 

and 

CORK  (A,  A  ,)  =  (II.D.2.4) 

n  n-i  k 

In  fact  An  is  just  a  GLAR(l)  process,  so  the  result  is  not 
surprising.  Using  II. 0.2. 2,  II.D.2.3,  and  an  induction  argu¬ 
ment  leads  to  the  general  m-step  correlation  formula 

CORR ( A  , A  )  =  n>m>  0.  (II.D.2.5) 

XI  XI  “Xu 

The  correlation  structure  between  {A  }  and  (G  }  can 

n  n 

be  derived  in  a  similar  fashion.  However,  it  is  more  direct 

to  note  that  since  the  {An>  sequence  is  the  same  as  the  GLAR(l) 

process  of  Section  II. B,  the  {An> ,  {Gn}  correlation  structure 

will  be  the  same  as  that  derived  for  the  {X  },  { G  }  sequences 

n  n 

in  II. B. 2.  Hence, 

C0RR(An,Gn_m)  =  (|)  (^)in,  n  >_  m  >  0.  (II.D.2.6) 

Of  course,  if  j  >  i;  then  C0RR(A^,G^)  =  0 

Now  the  serial  correlation  for  the  (Xn)  sequence  can 
be  found.  From  II. D. 2.1 
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Using  the  stationarity  of  the  (An)  sequence,  the  iid  nature 

and  independence  of  {B  },  {C  } ,  {G  } ,  the  fact  that  A.  is 

n  n  n  1 

independent  of  Gj  when  j  >  i,  and  the  fact  that  (An>  is  inde¬ 
pendent  of  { }  by  construction  and  taking  expectations,  we 
have 


E  (XX  ,)  =  E(B  )E(B  ,)E(A  ,A  ,)+E(Bn)E(C_  ,)E(A  ,C  ) 

n  n— l  n  n-l  n— 1  n-Z  n  n-l  n-1  n-l 


+  E  (B  ,  )  E  (C  )  E  (A  _ )  E  ( G  )  +E  ( C  )  E  ( C  ,)E(G„)E(Gri 

n-l  n  n-Z  n  n  n— 1  n  n-l 


=  [E  (B  )  ]  £  (A  -i  A  ,)+E(Bn)E(C_  ,  )  E  (A_  ,  G„  ,) 
n  n-l  n-z  n  n-l  n-l  n-l 


+  E(Bn_1)E(Cn)E(An_2)E(Gn)  +  [E(Cn)  ]  2  [E  (Gn)  ]  2  .  m 


(II .D.2. 7) 


From  II. D. 2.1 


=  [E (B  ) E (A  i ) +E(C  ) E (G  ) ]  x 
n  n-i  n  n 


[E  (B  .)E(A  ,)+E(C  .  )  E  (G  ,)], 

n— i  n-z  n— i  n— i 


i 


E (X  ) E (X  .) 
n  n— l 


E (X  ) E (X  , )  =  [E(B  )]2E(A  ,)E(A  -) +2E (B  ) E (C  ) E ( A  ) E (G  ) 

n  n-i  n  n~l  n-Z  n  n  n  n 

+  [E(Cn) ] 2 [E(Gn) ]2.  ( II .D .2.8) 

Using  II.D.2.7  and  II.D.2.8  to  compute  COV(Xn,Xn_^)  yields 

COV  (X  ,  X  ,)  =  [E  (B  )  ]  2  [E  ( A  A  ~)-E(A„  .  )  E  (A  ,)  ] 

n  n-i  n  n— i  n-z  n— l  n— z 

+  [E  (B  )  E  (C  )  ]  [E  (A  .)-E(A  )E(G  )] 

n  n  n-l  n-1  n  n 

Since  (A  },  (G  },  and  {X  }  are  all  Gamma(k,l),  VAR(X  )  =  VAR (A  ) 
n  n  n  n  n 

=  VAR(G  ) .  Hence, 
n 

CORK ( X, X  .)  =  [ E(B  ) ] 2CORR(A  , A  , ) + [E (B  ) E (C  ) ] CORK (A  , G  ) 
n  n— l  n  n  n—x  n  n  n  n 

From  II.D.2.4  and  II.D.2.6  we  know  this  equals 

CORR(Xn,Xn_1)  =  (^)2(^)  +  (^)  (|)  (f)  (II.D.2.9) 

Using  II.D.2.5  and  II.D.2.6  in  an  induction  argument 
yields  the  general  m-step  correlation  of 


GQRRU^Xn.m) 


=  2(^£)  +  (^)  (?)  (?)  ] 


k  '  k  vk' 


Recognizing  the  expression  in  brackets  as  CORR(Xn,Xn_^)  and 
letting  o  equal  this  correlation  we  have 
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CORR  (X  ,  X  J  =  ,  n  >  m  >  1. 

Xi  n- m  j\ 


(II. D. 2. 10) 


Figure  II. D. 2.1  shows  the  possible  combinations  of 
one-  and  two-step  correlations  for  the  GLARMA (1,1)  model. 

This  concludes  the  development  of  the  correlation  structure 
of  the  GLARMA (1,1)  model. 

3.  GLAR(2) 

Jacobs  and  Lewis  [Ref.  12]  first  developed  the  follow¬ 
ing  mixture  scheme  for  generating  a  p--order  autoregressive 
processes.  We  now  adopt  that  scheme  for  generating  a  second- 
order  autoregressive  sequence  of  random  variables.  This  is 
the  special  case  of  GLARMA (p,q)  with  p  =  2  and  q  =  0.  As  such 
it  closely  resembles  the  GLAR(l)  process.  Let 


X  =  B  X  +  CG  , 
n  n  n-Tn  n  n 


(II. D. 3.1) 


where  {Xn,  n  =  -1,0,1,...}  is  assumed  to  be  a  second-order 
stationary  sequence  of  Gamma (k,l)  random  variables;  {B^, 
n  =  1,2,...}  is  an  iid  sequence  of  Beta(k-q,q)  random  varia¬ 
bles;  { C^}  is  an  iid  sequence  of  Beta(q,k-q)  random  variables; 
{Gn}  is  an  iid  sequence  of  Gamma (k,l)  random  variables;  {Bn} , 

(C  },  {G  }  are  independent;  also  T  is  iid  with  ?  (_T  =1)  = 
n  n  n  n 

1  ~ P (Tn  = ^ )  =  u •  The  Gamma (k,l)  and  Beta(m,n)  densities  are 
found  in  II. B. 1.2  and  II.B.1.3,  respectively. 

This  generation  scheme  works  even  though  Xn_^  and 
Xn_2  are  dependent  random  variables.  The  mixture  of  Xn-1  and 
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Xn_2  produced  by  II. D. 3.1  is  Gamma  distributed  and  indepen¬ 
dent  of  Gn.  Hence,  II. D. 3.1  is  simply  another  example  of  the 
random  sum  of  two  independent  Gamma  random  variables  producing 
another  Gamma  random  variable. 

Two  special  cases  of  II. D. 3.1  are  as  follows.  When 
a*l,  the  GLAR ( 2 )  process  reduces  to  the  GLAR(l) .  When  q  =  k, 
the  { sequence  is  iid. 

The  serial  correlation  of  the  { }  sequence  can  be 
calculated  in  the  usual  fashion,  assuming  stationarity  of 
the  process  we  have 


=  B  X  m  +  C  G 
n  n-T  n  n 
n 


k>0;  0  <  q  <_  k 


and 


XX  .  =  B  X  X  .  +  C  G  X  .  , 

n  n-1  n  n-T  n-1  n  n  n-1 

n 


Using  the  independence  of  {Cn}  and  {G^} ,  the  fact  that  X^ 
is  independent  of  C ^ ,  G ^  and  B ^  when  j  >  i  and  taking  expec¬ 
tations  we  have 


E(XnXn-l)  =  ctE(Bn)E(x^_i)  +  (l-a)E(Bn)E(Xn-1Xn_2)+E  (C^E(Gn)E(Xn,1) 


-  ■><*?>  <^±i)  +  (1-0  <l£^a>E(Xn_1xn_2)  +|. 


since  E(Xn)  =  E(Gn)  =  1  by  assumption.  Using  the  second-order 
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stationarity  of  (Xn) 


E(XnXn-l)(1  “  (  — 9t^(-rq'-) 


_  a  (k-g)  (k+1)  +  g 
k2  k 


Upon  simplification 


E(Wi> 


_  gk  +gk-gkq-gq+kq 
k  (q+gk-gq) 


Hence, 


COV(X„,Xr,  ,) 
n  n— l 


_  i  i.\  r  ^  i 

k^  q+g  (k-q) 


and 


CO“>«„-Vl)  -  q+a (k-q)  (II.D.3.2) 

If  a  =  If  this  equals  1  -  if  k  =  q  then  it  is  zero.  Since 
q  >  k  it  is  clearly  non-negative. 

The  lag  two  correlation  can  be  calculated  in  a 
similar  fashion.  We  have 


Xn  "  BnXn-Tn+  CnV 


X  X„  -  =  B  X  _  X  _  +  C  G  X  _ 

n  n-2  n  n-T  n- 2  n  n  n-2 


and 
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E(Xnx„-2> 


-  a(!^)E(X„.1Xn.2)  +  (1-x)(^)E(xJ_2)  +3. 


Using  the  second-order  stationarity  of  the  { Xn>  sequence  we 
can  write 


E  (X  ) £  (X  -) 
n  n— 2 


=  a  (—3) E  (X  ^EtX^  ,)  +  (l-a)  (—3.)  [E(Xn  -)  ] 
q  n-i  n-2  K  n— 2 


so  that 


+  kIE(xn)l2' 


C0V(VXn-2> 


+  (l-o)  (^)  [E(XE_2)-{E(Xn_2)}2] 


+  I  *  f‘E(Xn)l2 


*  <*?>  1E  (xn-lXn-2>  -E  IXn-l)E (xn-2> 1 


+  (— a--^— [E(X^_2)  -  (E  cxn_2> } 2  ] 


“lT2>C0vlXn-l'X„-2)  +  ^^VAIUX^) 


Hence , 


CORK (X  , X  -) 
n  n—  2 


*  a(!Sia)COFR(Xn.1,XI1.2)  +  li^IOSral  , 
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CORR(Xn,Xn.2)  =  (ffi  (II.D.3.3) 

The  solution  procedure  for  II.D.3.3,  if  followed  for 
Xn-m'  w^-11  Produce  the  general  recursion  equation  (Yule-Walker) 
that  can  be  used  along  with  II.d.3.2  and  II.D.3.3  to  compute 
the  m-step  correlation.  The  formula  thus  produced  is 

CORR(Xn,Xn.m)  -  <*£2)  (aCORRlVVl-m’ 

+  (l-a)CORR(X  ,X  JJ,  (II.E.3.4) 

n  n+r— m 

n  >  m  >  1 . 


As  mentioned  in  previous  sections  the  (X  , G  )  pair 

n  n 

can  be  considered  to  be  a  correlated,  bivariate  pair  of 
Gamma (k,l)  random  variables,  Therefore,  we  proceed  to 
derive  the  correlation  structure  between  these  two  seauences . 
From  I I . D . 3 . 1  we  have 


Xn  =  BnXn-Tn+  CnGn 


and 


XG  =  BX  _  G  +  C  G  . 
n  n  n  n-TR  n  n  n 


Recalling  the  independence  of  { }  and  {Gn>  and 
when  j  >  i,  taking  expectations  yields 


X .  and  G . 
i  3 
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E(W 


<*£*>  +  (f) 


Hence, 


COV(Xn,Gn)  =  ^ 

and 


C0E*(W  *  ?• 


(II 


Continuing  this  process  we  have 


n 


B  X  _  +  C  G 
n  n-T  n  n 


and 


XG  , 
n  n-i 


=  Vn-TGn-l 


CG  G  1  * 
n  n  n-i 


Thus 


E  (X  G  -) 
n  n-1 


a(^3)E(Xn_1Gn_1)  +  (1-a)  (^2) 


+ 


2 

k 


Hence, 
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cov(VGn-l>  *  “(!S£2>cov(xn-l’Gn-l) 


And 


CORR(Xn,Gn_1)  =  a(^)(2). 


One  further  step  in  this  process  using  an  arbitrary 
m-step  lag  produces  the  general  recursion  formula 


CORR (X  , G  m) 
n  n-m 


(!s=2)  [aCORK(Xn,Gn+1.m) 


+  (l-a)CORR(X  ,G_ , -  ) ] ,  n  >  m  >  2  (II.D.3.7) 

n  n+  z— m  —  — 


Figure  II. D. 3.1  displays  a  plot  of  the  possible  com¬ 
binations  of  CORR ( X  . X„  .)  and  CORR(X  ,X  -) .  Note  that  when 

n  n— l  n  n— 2 

a  =  1,  the  GLAR(2)  process  reduces  to  the  GLAR(l)  and 

CORR(X  , X  -)  =  (^2)  2  which  defines  the  lower  boundary  of  the 

plot.  When  a  =  0,  CORR(X  .X  .)  =  0  and  CORR(X  ,X^  _)  = 

n  n-l  n  n-z  k 

which  goes  from  zero  to  one.  In  the  interior  of  the  graph, 
when  a  does  not  assume  an  extreme  value,  C0RR(Xn,Xn_2)  does 
not  reach  a  value  of  one.  This  is  demonstrated  by  the 
following  calculation.  From  II.D.3.3 


C0RRlxn'xn-2>  *  (Ta> 


If  this  correlation  is  to  equal  one,  then 
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Figure  II 


q+ak-2gq  _  k 
q+g (k-q)  -  k-q 


which  quickly  reduces  to 


q (-2gk-q+2kq)  =  0 


Since  0  <  q,  this  requires 


-2ak-q+2ctq  =  0 


(11.0.3.8 


Therefore, 


n  ~  _ —2 _ 

a  2 (q-k) ' 

Since  we  know  that  0  <  q  _<  k,  a  must  be  less  than  zero. 

However,  a  is  a  probability  and,  hence,  is  non-negative. 
Therefore,  the  original  requirement  in  II.D.3.8  cannot  be 
satisfied.  Hence,  C0RR(Xn,Xn_2)  cannot  equal  one. 

4 .  BGAR(l),  Bivariate  Model 

To  this  point  the  only  examples  of  bivariate  Gamma 
processes  presented  were  those  in  which  the  innovation  se¬ 
quence  was  one  part  of  the  bivariate  process  with  the  gener¬ 
ated  sequence  the  other  part.  The  simple  random,  linear 
structure  of  the  GLAR(l)  process  makes  it  easily  extendable 
to  a  variety  of  bivariate  models.  We  address  only  the 
simplest.  Consider  the  following  pair  of  random  variables, 
both  of  which  are  formed  from  the  same  innovation  process,  (Gn) 
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X„  =  B  X  ,  +  C  G  , 
n  n  n-1  n  n 


(II. D. 4.1) 


Y  =  D  Y  .  +  F  G  , 
n  n  n-1  n  n 


(II.D.4.2) 


where  {Xn,  n  =  0,1,...}  is  a  second-order  stationary  sequence 
of  Gamma (k,l)  random  variables;  {Yn,  n  =  0,1,...}  is  a 
second-order  stationary  sequence  of  Gamma (k,l)  random  varia¬ 
bles;  (Bn,  n  =  1,2,...}  is  an  iid  sequence  of  Beta(k-q,q) 
random  variables;  {Cn,  n  =  1,2,...}  is  an  iid  sequence  of 
Beta (q, k-q)  random  variables;  {Dn,  n  =  1,2,...}  is  an  iid 
sequence  of  Beta(k-r,r)  random  variables;  (Fn,  n  =  1,2,...} 
is  an  iid  sequence  of  Beta(r,k-r)  random  variable;  (Gn, 
n  =  1,2,...}  is  an  iid  sequence  of  Gamma (k,l)  random 

variables;  { B  } ,  (C  },  (D_  } ,  {F  },  and  (G  }  are  independent; 
n  n  n  n  n  ^ 

0  <  r  £  k;  0<q_<k. 

This  is  a  special  case  of  a  general  situation.  In 

a  more  general  case  the  {Xn}  and  {Yn}  sequences  could  have 

separate,  correlated  innovation  sequences  instead  of  sharing 

a  single  {G  }  sequence.  In  addition,  the  {B  }  and  {D  } 
n  n  n 

sequences  and  the  {Cn}  and  {F^}  sequences  could  be  correlated. 

Before  examining  the  correlation  between  (Xn}  and 
{Yn},  it  will  be  necessary  to  address  the  correlation  of  each 
of  these  sequences  with  the  {Gn}  sequence.  This  is  most 
easily  handled  by  recognizing  that  the  relationship  of  the 
{XR}  and  (Yn}  sequences  to  the  {Gn}  sequence  is  exactly  the 


same  as  the  (An}  sequence  to  the  CGn}  sequence  in  II.D.2.2. 
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Hence,  the  correlation  structure  will  be  the  same.  There¬ 
fore,  if  we  let  p  =  CORK (X  , X  . )  and  pv  =  CORR(Y  ,Y  .), 

a  a  n-x  x  n  a~  x 

n  =  1,2,...,  these  correlation  structures  can  be  written 
without  further  analysis  as 


CORR(X„.Gn)  =|=1-  CORR(Xn,Xn.1)  -  1  -  P*; 

C08R(VGn-l>  *  =  i1'  0X>  PX! 

C0RE(VGn-m>  =  (|)(^)“  =  (l-px)p”  n  >  m  >  0. 


(II.D.4.3) 

(II.D.4.4) 

(II.D.4.5) 


and 


CORR(Yn,Gn)  = 


CORR  (Y  , G  ,) 
n  n-x 

“«n'eJ 


I  *  1  -  C0RR<  Wl>  "  1  -  V 


-  <!><*?>  -  'i-w 

,  r.  ,k-rx  m  _  „  v  ni  „  ^  m  ^  „ 

Py)p  yf  ^  ^  ^  0 


(II. D. 4 .6) 

(II.D.4.7) 

(II.D.4.8) 


The  assumption  is  that  the  bivariate  process  is  stationary, 
although  it  should  be  noted  that  starting  the  univariate 
processes  in  a  stationary  mode  does  not  make  the  bivariate 
process  stationary.  The  initial  pair  {X0,YQ}  must  have  the 
bivariate  Gamma  distribution  associated  with  the  stationary 
Markovian  process . 

Now  we  address  the  cross  correlation  between  {X  }  and 

n 

{Yn>.  Start  with  II. D. 4.1.  We  have 
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n 


=  B  X  ,  +  C  G 
n  n-1  n  n 


and 


XY  ,  =  BX  ,Y  ,+CGY  , 

n  n-1  n  n-1  n-1  n  n  n-1 


Using  the  independence  of  {Cn>  and  {Gn>  and  that  of  Y 
Gj  and  Cj  when  j  >  i  and  taking  expectations 


E(X  Y  )  =  E (B  ) E (X  i Y  .)  +  E (C  ) E ( G  ) E ( Y  , 

n  n-l  n  n-l  n-1  n  n  n— 1 


so  that 


E(X  Y  ,) 
n  n-1 


*  +  k- 


Now  starting  with  II.E.4.2,  we  have 


Y  =  D  Y  ,  +  F  G 
n  n  n-1  n  n 


and 


X  Y 
n  n 


=  D  Y  .X  +  F  G_X  , 
n  n-1  n  n  n  n 


Taking  expectations  as  above 


E(x  Y  }  =  E (D  ) E (X  Y  ,)  +  E (F  ) E ( G  X  ) . 

n  n  n  n  n- 1  n  n  n 


and 


)  , 


(II.D.4.9) 
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(II. D. 4. 10) 


Using  II.D.4.3  we  deduce  E(G  X  )  =  —  so  that 

n  n  ^ 


E,Vn'  *  +  ‘fxM3’ 

k 


Invoking  the  second-order  stationarity  of  the  (Xn,Yn)  se¬ 
quences,  using  II.D.4.9  and  II. D. 4. 10  and  letting  w  =  E (X^Y^) 
and  z  =  ECX^),  we  have  the  two  equations 


w 


<#> 


(II. D. 4 


z 


(II.  D. 4 


Using  II. D. 4. 12  to  substitute  for  z  in  II. E. 4. 11  yields 

»  -  <¥>««¥>- +  l> +  <!><#> 

k 

=  r>s2r.kg-lsq+cirlw  +  (k~r)g  (k2,t9L).r 

1  u2  J  .2  ,3 

k  k  k 

When  E (XnYn)  is  substituted  for  w  after  simplification,  this 
produces 


E(XnYn> 


2  2 
k  q-kqr+k  r+rq 

k(kq+kr-qr) 


Hence 


COV (X  , Y  ) 
n  n 


rq 

k (kq+kr-qr) 


(II. D. 4 


.11) 

.12) 


.13) 
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and 


CORR(Xn,Yn) 


kq+kr-qr 


,  k>0,  0<q<k,  0  <  r  <  k . 


(II. D. 4. 14) 


(1-p  ) (l-pY> 

_ X _ I  - 

(1-p  )+(l-p  ) p 


(1-PX) (l-py) 

(1-PY)+(1-PX) PY  * 


(II. D. 4. 15) 


This  latter  expression  follows  since  for  a  given  k  the  corre¬ 
lation  structure  is  parameterized  by  r  and  q  or  equivalently 
by  the  serial  correlations  px  and  pY  given  at  II.D.4.3  and 


II.D.4.6 


We  can  now  substitute  II. E. 4.13  into  D.E.4.12  to 


solve  for  E(X  Y  .) 

n  n-i 


E(X  Y 
n  n-i 


2  2 

,k-q.  .k  g-kqr+k  r+rq-. 

'  '  t  Ir  rirr.il*-. v- \  J 


k  (kq+kr-qr) 


3  3  2  2 

k  q+k  r+kqr-k  qr-q  r 


k  (kq+kr-qr) 


Hence , 


“'"Wl1  * 


kqr-q  r 

2 

k  (kq+kr-qr) 


CCRR(X  ,Y  . ) 
n  n-i 


( _ 2£ _ )  (]£i2) 

vkq+kr-qr;  v  k  ' 


CORR(X  Y  ) p 
n  n  X 


(II. D. 4 .16) 
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Continuing  in  this  vein  produces  the  general  formula 


CORS(Wm> 


,_qr_.  ,k-q,  m 
vkq+kr-qr'  v  k  ' 


=  CORR(Xn/Y^  p™,  n  >  m  >  0. 


(II. D. 


Solving  for  correlations  where  the  {Xn)  lags  the 
{ Yn }  sequence  is  similar  to  the  above  process,  but  somewhat 
abbreviated.  Starting  with  II.D.4.2, 


n 


=  D  Y  i  +  F  G  , 
n  n-1  n  n 


we  have 


Y  X  . 
n  n-1 


=  D  Y  .  X  ,  +  FGX 

n  n-1  n-l  n  n  n-1 


Taking  expectations  as  before  gives 


E<Vn-l> 


E ( D  ) D ( Y  ,Xn  .) 
n  n-l  n-l 


E<Fn)E(Gn> 


E(xn-1> 


Using  the  second-order  stationarity  of  (Xn)  and  {Yn}, 

E (Y  .X  . )  is  known  from  II. E. 4. 13,  so 
n-l  n-l 


E(YnX 


n-1 


2  2 

,k-rs  ,k  q-kqr+k  r+rq,  ,  r 
(~TT)  (  ~k  [kq+kr-q'r  ]  +  k 


k^q+k2r+kqr-k2qr-qr2 
k2 (kq+kr-qr) 


.17) 


Hence, 


COVlVYn-l> 


.  2 
kqr-qr 

(kq+kr-qr) 


and 

C0EE(Wl>  *  (gq+l-qr1  (!T£l  *  C01«(X„Yn)pY. 

Further  computations  of  this  nature  produce  the 
general  formula 

CORE(Yn,Xn.m)  -  (*?>" 

=  CORlUX^Y^  p”  n>m>0;k>0;0<q<k; 

0  <  r  5  k. 

To  examine  the  correlation  further,  note  that  if 
ox  =  PY  =  P>  then  II. D. 4. 15  yields 

CORR(Xn,Yn)  .  &JL. 

Thus  if  o  -  0  (i.e.,  the  {X  }  and  {Y  }  processes  are  iid 

n  n 

sequences) ,  this  correlation  is  one.  For  {X  }  and  {Y  }  to 

n  n 

be  iid,  we  must  have  X  =  G  =  Y  .  Therefore,  this  limiting 

n  n  n 

correlation  does  make  sense  and  suggests  that  perhaps  a  bi¬ 
variate  Gamma  should  be  used  as  the  innovation  process.  This 
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would  allow  for  separate  control  of  the  serial  correlation 

(auto-correlation)  and  cross-correlation  of  the  (X  }  and  {Y  } 

n  n 

sequences . 

If  p  +  1  (i.e.,  Xn  ~  Xn-1'  Yn  ~  Yn-1^ '  the  e^fect  of 
the  innovation  sequence  is  slight  and  the  cross-correlation 
between  {Xn}  and  (Yn>  goes  to  zero.  In  a  more  complicated 
model  than  we  have  addressed  here  the  cross-correlation  may 
be  controlled  by  imposing  a  correlation  on  the  {Bn>  and  {Dn> 
sequences . 

Cross-coupled  processes,  as  discussed  in  Gaver  and 
Lewis  [Ref.  2] ,  are  possible.  These  can  be  used  to  create 
negative  serial  correlations  in  the  { }  and  (Y  }  processes. 

E.  NUMERICAL  CONVERGENCE  OF  THE  MAXIMUM  LIKELIHOOD  COMPUTER 
PROGRAM  AND  SIMULATION  STUDY  OF  PROPERTIES  OF  ESTIMATORS 

The  program  described  in  Section  II. B. 4  for  computing  the 

conditional  density  function  in  the  GLAR1  process  was  used 

in  two  fashions.  First,  it  was  tested  by  using  computer 

generated  data  from  a  GLAR(l)  process  with  known  parameter 

values  k,  q,  and  u-  Simultaneously  a  simulation  of  properties 

of  m.l.e  's  k  and  q  for  k  and  q  was  conducted.  Second,  it  was 

used  to  estimate  the  parameters  in  the  GLAR(l)  model  for  real 

wind  speed  data.  Only  the  first  use  is  covered  in  this  section. 

The  second  use  is  addressed  in  Chapter  IV,  Preliminary  Data 

Analysis . 

Four  aspects  of  the  program  were  addressed  in  the  use  of 

the  program  with  simulated  data.  These  were:  sensitivity 

of  the  maximum  seeking  method  to  start  point,  the  size  of  the 
♦maximum  likelihood  estimate 
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standard  deviation  of  the  m.l.e  and  moment  estimates  of  k  and 
q  produced,  and  the  degree  of  bias,  if  any,  in  the  estimates. 

In  addition,  normality  of  the  distributions  of  the  estimates 
was  investigated  using  normal  plots  and  Shapiro-Wilks  tests. 

Because  of  the  large  computation  time  involved  in  obtain¬ 
ing  a  m.l.e* the  simulation  study  was  small.  Three  types  of 
data  generated  from  GLAR(l) processes  were  used  to  exercise  the 
program.  Each  type  of  data  consisted  of  ten  independent  sets 
(replications)  of  1000  data  points  each.  The  k  and  q  values  and 
the  correlation  were  varied  from  one  type  of  data  to  another. 
Thus  the  first  type  of  data  was  generated  with  a  k  value  of 
4.0  and  a  q  value  of  1.0.  These  parameter  values  produce  a 
correlation  of  0.75  (see  equation  II. B. 2.1).  The  second  type 
of  data  varied  the  correlation,  but  retained  the  same  k  value. 

A  k  of  4.0  and  a  q  of  3.0  produce  a  correlation  of  0.25.  These 
values  were  used  to  produce  data  sets  of  the  second  type. 

The  third  type  of  data  returned  to  a  high  correlation,  but 
used  a  small  k  value.  The  parameter  values  used  to  generate 
this  data  type  were  a  k  of  0.75  and  a  q  of  0.1875.  These 
values  also  produce  a  correlation  of  0.75.  Table  II.E.l 
summarizes  these  cases.  In  all  cases,  y  =  1. 

TABLE  II.E.l 


CASE 

k 

q 

P 

1 

4.0 

1.0 

0.75 

2 

4.0 

3.0 

0.25 

3 

0.75 

0.1875 

0.75 

*maximum  likelihood  estimates 


The  Gamma  variates  were  generated  by  the  program  LLRANDOMII 


[Ref.  16]  and  all  runs  were  performed  on  the  NPS  IBM/3033 
computer. 

To  test  the  sensitivity  of  the  maximum  likelihood  computer 

program  to  the  starting  point  of  the  search  for  a  maximum, 

each  set  of  data  was  used  in  two  runs  of  the  program.  The 

first  run  used  the  actual  parameter  values  k  and  q  as  a  start 

point.  The  second  run  used  the  moment  estimates  k  and  q  of 

the  parameter  values  (see  equations  II. B. 4. 15  and  II. B. 4. 16)  as 

*  ~ 

a  start  point.  The  resulting  m.l.e.  parameter  estimates  k  and 

/s. 

q  were  recorded. 

The  first  case  had  k  =  4.0  and  q  =  1.0.  The  results  of 
the  computer  runs  are  presented  in  Table  II. E. 2  for  data  of 
the  first  type.  The  last  column  presents  the  two-dimensional 
distance  in  the  (k,q)  plane  between  the  estimates  produced  by 
the  two  different  start  points  for  each  data  set.  Although 
the  values  do  not  differ  widely,  the  relatively  large  differ¬ 
ences  in  some  cases  indicate  that  the  likelihood  function  is 
relatively  flat  near  the  maximum. 

However,  there  is  another  factor  which  may  be  contributing 
to  differences  in  final  parameter  estimates.  The  calculation 
of  the  likelihood  function  for  1000  data  points  requires  the 
numerical  evaluation  of  999  integrals  (see  equations  II.B.4.2 
and  II.B.4.3).  The  IMSL  routine  D CAD RE  was  used  for  this 
evaluation.  Two  of  the  parameters  in  the  call  to  DCADRE  are 


* 


maximum  likelihood  estimates. 
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TABLE  II. E. 2 


Results  of  Search  for  Maximum  Likelihood  Estimates 
in  a  GLAR(l)  Model;  k  =  4.0,  q  =  1.0 


Data 

Set 

Starting 

k 

Value 

q 

Run  Time 
(in  Min.) 

Nurrber  of 
Iterations 

Ending  Value  Iteximum 
k  a  Likelihood 

0 

4.000 

1.000 

20 

6 

3.668 

0.882 

-214.479 

0 

2.884 

0.538 

129 

17 

3.624 

0.886 

-214.488 

1 

4.000 

1.000 

19 

2 

4.004 

0.979 

-204.819 

1 

4.235 

1.101 

46 

5 

4.014 

0.983 

-204.820 

2 

4.000 

1.000 

104 

14 

3.451 

0.869 

-260.332 

2 

3.360 

0.805 

112 

14 

3.440 

0.867 

-260.333 

3 

4.000 

1.000 

16 

2 

3.977 

1.051 

-206.416 

3 

3.767 

0.947 

71 

15 

3.955 

1.041 

-206.418 

4 

4.000 

1.000 

52 

5 

4.246 

1.232 

-304.365 

4 

4.423 

1.233 

45 

6 

4.254 

1.235 

-304.366 

5 

4.000 

1.000 

47 

3 

4.061 

0.998 

-211.074 

5 

4.647 

1.187 

27 

4 

4.122 

1.028 

-211.060 

6 

4.000 

1.000 

61 

4 

3.593 

0.891 

-229.087 

6 

3.387 

0.782 

55 

9 

3.565 

0.883 

-229.092 

7 

4.000 

1.000 

37 

5 

4.041 

0.973 

-241.637 

7 

4.421 

1.069 

82 

7 

4.051 

0.974 

-241.636 

8 

4.000 

1.000 

37 

6 

4.024 

0.999 

-240.432 

8 

3.927 

0,856 

24 

8 

4.106 

1.033 

-240.424 

9 

4.000 

1.000 

51 

3 

4.109 

1.054 

-255.245 

9 

4.398 

1.112 

38 

5 

4.174 

1.081 

-255.229 

Value  of 
Difference 
(10“3) 

46 

10 

11 

? 

24 

8 

68 


29 


10 


88 

70 


the  relative  and  absolute  errors  allowed  for  this  calculation. 

Practical  considerations  of  computer  run  time  dictated  rela- 

-4 

tively  large  values  of  10  for  each  of  these  parameters. 

This  error  when  applied  999  times  in  the  calculation  of  the 
likelihood  function  may  have  lead  to  the  differences  in  m.l.e. 
parameter  estimates.  As  a  test  of  this  hypothesis  the  data 
set  which  produced  the  largest  distance  between  the  pairs  of 
estimates  (data  set  8)  was  rerun  with  DCADRE  error  parameters 
set  at  10  The  run  which  used  the  actual  parameter  values 

as  a  start  point  ran  in  171  minutes  and  produced  estimates  of: 

A  A 

k  =  4.102,  q  =  1.031.  The  run  which  used  the  moment  estimates 
as  a  start  point  was  terminated  after  404  minutes  CPU  time.  At 

A  /V 

that  point  it  had  estimates  of:  k  =  4.088,  q  =  1.025.  The 

_3 

distance  of  15  x  10  represents  a  significant  reduction  in 
the  previous  distance  of  88  *  10-3.  It  seems  from  this  exam¬ 
ple  that  the  program  can  be  made  less  sensitive  to  the  starting 
point  by  increasing  the  accuracy  with  which  DCADRE  computes 
the  integrals  in  the  likelihood  function.  Of  course,  a  con¬ 
siderable  increase  in  computational  cost  is  incurred.  This 
cost  is  not  practical  in  these  simulations  or  necessary  since 
only  a  rough  idea  of  the  properties  of  the  estimates  was  sought. 
The  results  of  the  runs  for  data  type  one  are  also  presented 

^  A  A 

in  Figure  II.E.l.  First,  each  pair  of  estimates,  (k,q)  and  (k,q) 
is  plotted  in  the  k,q  plane.  Then  each  point  is  projected 
along  each  axis  to  more  conveniently  reflect  the  marginal 
variation  in  the  estimates  for  each  parameter. 
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MOMENT  ESTIMATES  (O) 


3.50  4.00  4.50 

VALUE 


When  k  =  4.0  and  q  =  1.0,  the  method  of  moments  produces 
estimates  for  k  with  a  sample  mean  of  3.94  and  a  sample 
standard  deviation  of  0.55.  The  statistics  for  the  estimates 
of  q  are  a  sample  mean  of  0.96  and  a  sample  standard  deviation 
of  0.21.  The  maximum  likelihood  method  produces  estimates  of 
k  with  a  sample  mean  of  3.93  and  a  sample  standard  deviation 
of  0.27.  The  values  for  q  estimates  are  a  sample  mean  of  1.00 
and  a  sample  standard  deviation  of  0.11.  Although  the  method 
of  moments  and  maximum  likelihood  method  do  not  produce  signi¬ 
ficantly  different  values  for  the  mean  of  the  estimates  of  the 
parameters,  the  lower  estimated  standard  deviation  for  the 
likelihood  method  makes  this  technique  more  desirable  from 
the  standpoint  of  precision.  No  bias  was  evident  in  either 
estimation  technique  with  the  precision  attained  in  the 
simulations . 

The  second  case  was  the  low  correlation  case  with  k  =  4.0 
and  q  =  3.0.  Here  the  distinction  between  the  two  estimation 
procedures  is  not  as  clear  (Table  II. E. 3  and  Figure  II. E. 2). 

The  method  of  moments  produced  estimates  for  k  with  a  sample 
mean  of  3.99  and  sample  standard  deviation  of  0.17.  The  esti¬ 
mates  for  q  had  a  sample  mean  of  2.9T  and  sample  standard 
deviation  of  0.18.  The  maximum  likelihood  method  produced 
estimates  of  k  which  had  sample  mean  4.00  and  sample  standard 
deviation  0.16.  The  q  estimates  had  a  sample  mean  of  2.99 
and  sample  standard  deviation  of  0.17.  It  is  clear  that  neither 
method  of  parameter  estimation  enjoys  a  distinct  advantage 
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TABLE  II. E.  3 


Results  of  Search  for  Maximum  Likelihood  Estimates 
in  a  GLAR(l)  Model;  k  =  4.0,  q  =  3.0 


Value  of 


Data 

Set 

Starting 

k 

Value 

q 

Run  Time 
(in  Min.) 

Number  of 
Iterations 

Ending 

k 

Value 

q 

Maximum 

Likelihood 

Difference 

do-3) 

0 

4.000 

3.000 

38 

6 

4.078 

3.187 

-620.637 

8 

0 

4.019 

3.271 

25 

5 

4.085 

3.192 

-620.637 

1 

4.000 

3.000 

23 

4 

3.780 

2.934 

-655.488 

5 

1 

3.718 

2.843 

26 

4 

3.775 

2.931 

-655.489 

2 

4.000 

3.000 

40 

3 

3.789 

2.765 

-663.886 

3 

2 

3.912 

2.862 

27 

5 

3.789 

2.762 

-663.886 

3 

4.000 

3.000 

44 

7 

4.321 

3.132 

-585.776 

5 

3 

4.000 

3.116 

50 

4 

4.316 

3.130 

-585.774 

4 

4.000 

3.000 

44 

4 

4.166 

3.403 

-604.719 

22 

4 

3.958 

3.174 

68 

6 

4.152 

3.385 

-600.715 

5 

4.000 

3.000 

27 

6 

3.890 

2.904 

-638.667 

5 

4.001 

3.048 

13 

4 

3.885 

2.902 

-638.666 

5 

6 

4.000 

3.000 

43 

3 

4.051 

2.873 

-599.581 

6 

4.221 

2.902 

49 

5 

4.049 

2.878 

-599.581 

5 

7 

4.000 

3.000 

23 

2 

3.954 

2.934 

-627.062 

11 

7 

4.079 

3.074 

16 

7 

3.963 

2.941 

-627.058 

8 

4.000 

3.000 

21 

4 

3.917 

2.898 

-637.234 

11 

8 

3.708 

2.633 

32 

6 

3.908 

2.891 

-637.233 

9 

4.000 

3.000 

15 

2 

4.111 

2.924 

-590.563 

14 

9 

4.093 

2.906 

43 

2 

4.121 

2.934 

-590.562 

E.  (+)  AND  MOMENT  ESTIMATES  (©) 
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over  the  other  with  respect  to  precision  or  accuracy.  How¬ 
ever,  the  method  of  moments  is  considerably  cheaper  with 
respect  to  computation  time  required.  This  is  consistent 
with  known  results  that  for  i.i.d.  Gamma  data  (k  =  q) ,  the 
moment  estimate  for  k  is  quite  efficient  when  compared  to 
the  maximum  likelihood  estimator. 

Third  case  (k  =  0.75,  q  =  0.1825).  The  third  type  of  data 
was  a  high  correlation  case  with  a  low  k  value.  Specifically 
k  =  0.75,  q  =  0.1825,  and  the  correlation  was  0.75.  As  can 
be  seen  in  Figure  II. E. 3  and  Table  II.E.4,  both  the  methods  of 
parameter  estimation  considerably  overestimated  the  parameter 
values,  indicating  considerable  bias  in  the  procedures.  The 
method  of  moments  produced  estimates  for  k  with  sample  mean 
of  0.8061  and  sample  standard  deviation  of  0.067.  Those  for 
q  had  a  sample  mean  of  0.232  and  sample  standard  deviation  of 
0.026.  The  likelihood  method  produced  estimates  for  k  with  a 
sample  mean  of  0.852  and  sample  standard  deviation  0.039.  The 
corresponding  statistics  for  q  estimates  are  0.224  and  0.004. 
As  was  true  in  the  other  high  correlation  case,  the  standard 
deviations  of  the  maximum  likelihood  estimates  are  consider¬ 
ably  smaller  than  those  of  the  moment  estimates.  However, 
since  the  evidence  is  that  the  estimates  are  highly  biased, 
the  advantage  of  this  smaller  standard  deviation  is  not  clear 
unless  additional  data  would  serve  to  reduce  the  apparent 
bias.  The  detailed  results  for  this  data  type  are  presented 
in  Table  II.E.4  and  Figure  II. E. 3.  It  would  be  of  interest 


TABLE  II. E. 4 


Results  of  Search  for  Maximum  Likelihood  Estimates 
in  a  GLAR(l)  Model;  k  =  0.75,  g  =  0.1875 


Value  of 


Data 

Set 

Starting 

k 

Value 

q 

Run  Time 
(in  Min.) 

Number  of 
Iterations 

Ending 

k 

Value 

q 

Maximum 

Likelihood 

Difference 

(10-3) 

0 

0.7500 

0.1825 

85 

5 

0.861 

0.221 

325.685 

o 

0 

0.9795 

0.2489 

136 

4 

0.863 

0.222 

325.685 

z 

1 

0.7500 

0.1825 

59 

4 

0.832 

0.235 

866.216 

n 

1 

0.8026 

0.2351 

19 

1 

0.832 

0.235 

866.214 

(J 

2 

0.7500 

0.1825 

72 

4 

0.806 

0.224 

439.198 

1 

2 

0.7621 

0.2147 

66 

4 

0.807 

0.224 

439.198 

X 

3 

0.7500 

0.1825 

141 

2 

0.866 

0.222 

1131.684 

3 

0.7514 

0.2213 

121 

1 

0.866 

0.222 

1131.671 

0 

4 

0.7500 

0.1825 

59 

4 

0.807 

0.223 

2333.130 

4 

0.7749 

0.2420 

54 

5 

0.807 

0.222 

2333.131 

1 

5 

0.7500 

0.1825 

59 

4 

0.910 

0.217 

883.801 

5 

0.8336 

0.2320 

210 

4 

0.910 

0.217 

883.801 

0 

6 

0.7500 

0.1825 

129 

5 

0.885 

0.222 

287.798 

6 

0.7308 

0.1790 

88 

4 

0.884 

0.222 

287.799 

1 

7 

0.2500 

0.1825 

71 

5 

0.795 

0.223 

977.265 

7 

0.7759 

0.2241 

117 

8 

0.795 

0.223 

977.265 

0 

8 

0.7500 

0.1825 

66 

6 

0.906 

0.227 

1018.470 

8 

0.8092 

0.2863 

136 

6 

0.906 

0.227 

1018.472 

0 

9 

0.7500 

0.1825 

57 

4 

0.818 

0.212 

410.842 

36 

9 

0.8335 

0.2347 

85 

3 

0.852 

0.225 

411.703 
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to  see  if  a  technique  such  as  (2-fold)  jacknifing,  such  as 
that  applied  by  Quenouille  [Ref.  17]  to  correlation  estimates, 
would  help  here. 


A  much  larger  study  would  be  needed  to  come  to  definite 
conclusions  about  the  efficiency  of  maximum  likelihood  esti¬ 
mation  in  this  model.  However,  as  in  the  i.i.d.  case  for 
Gamma  variates,  m.l.e.'s  are  likely  to  be  considerably  better 
than  moment  estimations  for  values  of  k  less  than  one. 

Note  too  that  the  maximum  likelihood  estimation  did  not 
include  the  mean  value  parameter.  This  could  be  done  or  the 
mean  could  be  estimated  from  the  sample  mean  X.  The  inflation 
of  variation  of  X  due  to  the  correlation  is  known  to  be 
(asymptotically) 


1  +  2  l  pk 

k=l 


1  +  P 
1  -  P 


Thus  for  p  =  0.75,  the  effective  sample  size  for  estimating 
y  in  a  sample  size  of  size  n  is  n (1-p ) / (1+p) .  For  p  =  0.75 
this  is  n/7. 

The  normality  of  the  estimates  was  investigated  with  normal 
plots  and  Shapiro-Wilk  tests  for  normality.  A  summary  of 
results  if  given  in  Table  II. E. 5.  The  normality  hypothesis 
is  accepted  at  a  0.95  level  if  the  Shapiro-Wilk  statistic  W 
is  higher  than  0.842,  at  a  0.99  level  level  if  it  is  higher 
than  0.781.  No  strong  indication  of  non-normality  is  indicated 


Summary  of  Simulation  Results 
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III.  MOVING  AVERAGE  MODELS 

A.  INTRODUCTION 

Although  several  researchers  have  proposed  models  for 
correlated,  marginally  exponential  random  sequences  [Refs.  18, 
19,201,  Gaver  and  Lewis  [Ref.  2]  produced  the  first  analytically 
and  computationally  tractable  model  for  the  generation  of 
correlated,  marginally  Exponential  random  sequences.  They 
showed  that  in  the  usual  linear,  additive,  first-order,  auto¬ 
regressive  equation 

X  =  0Xn..  +  E  (III.A.l) 

n  n+1  n 

where  {Xn>  n  =  0,1,2,...}  is  a  second-order  stationary, 
marginally  Exponentially  distributed  sequence  of  random  varia¬ 
bles;  {En,  n  =  1,2,...}  is  an  independent,  identically  dis¬ 
tributed  (iid)  sequence  of  innovative  random  variables; 

0  <_  3  <  1,  the  distribution  of  the  {E^}  which  produces  the 
desired  marginal  distribution  for  {Xn}  is 

(0  with  probability  8, 

E  =  '  (III. A. 2) 

n  |  En  with  probability  1-6, 

where  {En,  n  *  1,2,...}  is  an  iid  sequence  of  Exponentially 
distributed  random  variables  with  the  same  parameter.  A,  as 
the  {Xn}  sequence.  Equation  III.A.l  can  now  be  written  as 


99 


with  probability  0 


(III. A. 3) 


6X 


n-1 


"  'I 


0Xn-1  +  En  with  probability  1-0 


If  {In>  n  =  1,2,...}  is  an  iid  binary  sequence  independent  of 
{Xn}  and  {En>  such  that  P(In  =  0)  =  1  -  P(In  =  1)  =  0,  then 
equation  III. A. 3  can  be  more  succinctly  written  as 


n 


=  0Xn  .  +  IE, 
n-l  n  n 


(III .A. 4) 


Xn  is  a  random  linear  combination  of  identically  distributed 

random  variables  in  the  sense  that  the  variable  En  actually 

enters  into  the  sum  only  when  the  random  variable  I  has  value 

one.  Since  for  a  given  {In}  sequence,  the  distribution  of  Xr 

depends  only  on  the  distribution  of  X  .  and  E  ,  X  will  be 

n-l  n  n 

Exponentially  distributed  whenever  both  X  .  and  E  are  inde- 

n- 1  n 

pendent  and  have  the  same  Exponential  distribution.  This 
understanding  allows  the  autoregressive  relation  in  III. A. 4 
to  be  transformed  into  a  moving  average  by  substituting  another 
innovative  random  variable  for  the  Xn_^  to  produce 


=  SE 


n 


InEn+l ' 


(III. A. 5) 


This  model  was  designated  the  EMA(l)  for  Exponential  moving 
average  of  order  one.  This  EMA(l)  model  is  one  dependent  in 
that  X^  and  xn+j  are  independent  for  j  ?  ±1.  Consequently, 
only  the  lag  one  correlation,  p1  =  =  CORR(Xn  ,X  ,  or  more 

completely  only  the  joint  distribution  of  Xn  and  Xn+^  need  be 
studied. 
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Lawrance  and  Lewis  [Ref.  5]  give  a  fairly  complete  des¬ 


cription  of  this  EMA(l)  model.  Of  particular  note  is  the 
relative  tractability  of  this  model  which  enabled  the  authors 
to  derive  correlations,  distributions  of  sums  of  Xn's, 
intensity  function,  spectrum  of  counts,  joint  density  of 
Xn  and  Xn+1,  conditional  expectations,  and  other  properties. 

The  existence  of  these  characteristics  is  beneficial  in  data 
analysis  and  is  a  primary  advantage  of  the  EMA(l)  over  pre¬ 
viously  suggested  models.  However,  the  EMA(l)  model  does 
possess  a  limitation.  The  range  of  possible  positive  correla¬ 
tions,  p^,  is  restricted  to  the  interval  from  zero  to  one- 
quarter.  Thus  for  a  given  correlation  between  zero  and 
one-quarter,  the  structure  of  Xn  and  Xn+^  and  the  sample  path 
behavior  of  the  sequence  are  determined. 

The  structure  of  III.A.l  is  that  of  a  special  random  linear 
combination  of  Exponential  random  variables  to  given  an  Exponen¬ 
tial  random  variable.  Other  such  random  linear  combinations 
are  now  known  and  for  the  first-order  autoregressive  case 
produce  dramatic  differences  in  sample  path  behavior  of  the 
sequence  {Xnl.  In  this  section  of  the  thesis  we  investigate 
these  random  linear  combinations  in  the  context  of  the  first- 
order  moving  average  structure. 

In  this  Chapter  we  examine  extensions  of  the  model  in  four 
ways : 

1 .  Negative  Correlation 

McKenzie  [Ref.  21]  has  suggested  a  scheme  for  inducing 
negative  correlation  in  the  EMA(l)  process  by  correlating  the 


I 


sequence  {I  } .  A  better  scheme,  in  that  it  produces  a  larger 
range  of  negative  correlations,  was  introduced  by  Lawrance 
and  Lewis  [Ref.  8]  for  the  extended  first-order  autoregressive 
model  NEAR { 1 ) .  This  scheme  involves  a  bivariate  error  se¬ 
quence  (E  , E^}  and  its  use  is  investigated  in  this  thesis  for 
the  moving  average  process . 

2.  New  Exponential  Movinq  Average  Model  of  Order  One, 

NEMA(l) 

It  will  be  shown  that  no  first-order  moving  average 
process  which  is  a  random  linear  combination  of  Exponential 
random  variables  can  have  correlation  greater  than  one-quarter. 
Thus  the  differences  in  the  processes  for  given  correlation  is 
investigated  in  terms  of  the  joint  structure  of  Xn  and  Xn+^. 

The  NEMA(l)  process  obtained  by  using  the  NEAR ( 1 )  structure 
[Ref.  8]  in  a  moving  average  context  is  analytically  tractable 
and  quantities  such  as  the  joint  Laplace-Stilt jes  transform 
of  Xn  and  Xfl+1,  the  spectrum  of  counts,  P(Xnfl  >  X^  ,  and  condi¬ 
tional  expectations  are  obtained.  It  also  combines  the  forward 
and  backward  EMA(l)  models  as  extreme  cases  and  is  thus  a 
natural  extension  of  the  EMA(l)  model. 

3 .  The  Moving  Minimum  Model 

A  non-linear  combination  of  Exponentials  involving 
minima  has  been  applied  by  Tavares  [Ref. 22]  to  obtain  a  first- 
order  autoregressive  process  which  is  intimately  connected 
[Ref . 23]  with  the  EAR(l)  process  of  Gaver  and  Lewis  [Ref.  24]. 
This  structure  is  applied  to  the  moving  average  process.  It 
is  shown  that  this  process  extends  the  range  of  attainable 
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correlations  in  first-order  "moving  average"  process  beyond 
that  obtainable  by  random  linear  combinations  of  Exponentials , 
However,  the  process  is  slightly  degenerate  in  that  there  is 
a  positive  probability  that  successive  points  will  lie  on  the 


line  =  bXn  for  positive  correlations  and  on  the  curve 

Xn  — ^n^l 

e  +  e  =1  for  negative  correlations.  The  price  paid 


for  the  extended  range  of  correlations  is  a  limited  analytical 
tractability  as  compared  to  the  EMA(l)  or  NEMA(l)  processes. 

The  moving  minimum  process  is  investigated  in  terms  of  the 
joint  structure  of  Xn  and  Although  the  joint  distribution 

can  be  derived,  the  functions  are  difficult  to  examine  in 
detail.  Therefore,  simple  characterizations  of  the  joint 
structure,  in  addition  to  the  correlation,  are  examined.  These 
include  the  P(Xn+1  > Xn) ,  a  crude  measure  of  the  tendency  of 
the  sequence  to  exhibit  runs  up  and  down,  and  conditional 
expectations,  E  (xn !  xn_]_  =  x)  and  E  (xn_1 !  xn  =  y)  - 
4 .  The  Beta-Exponential  Model 

Finally,  another  random  linear  combination  of  Exponen¬ 
tials  to  produce  correlated  Exponentials  is  examined.  Unlike 
the  previous  models,  the  coefficients  of  the  Exponential  random 
variables  are  themselves  continuous  random  variables.  This  in¬ 
creases  the  complexity  and  reduces  the  analytical  tractability 
of  this  model.  Simple  sample  path  characteristics  are  derived 
or  simulated.  These  are  special  cases  of  the  GLyiA(l)  from  Chapter  II 


B.  NEGATIVE  CORRELATION  IN  MOVING  AVERAGE  MODELS 

The  problem  of  non-negative  correlations  was  addressed  by 
McKinzie  [Ref.  21]  who  modified  the  form  of  the  EMA(l)  model  to  be: 


.03 


(III.B.l) 


n 


=  SE„  +  (1  - 


n 


Z  )  E  ,  , 
n  n- 1 


where  (En,  n  =  0,1,2,...}  is  an  iid  sequence  of  Exponential 
random  variables,  {Zn,  n  =  1,2,...}  is  a  sequence  of  binary 
random  variables  with  P(Zn=l)  =  1  -  P(Zn=0)  =  3,  and  Zn 
is  independent  of  all  En  and  all  past  Xn-  By  imposing  an 
MA (1)  correlation  on  the  sequence  {ZnJ,  McKenzie  was  able  to 
produce  a  negative  correlation  for  the  {Xn} .  However,  this 
negative  correlation  is  achieved  at  the  cost  of  reducing  the 
possible  range  of  positive  correlations  for  the  {Xn}.  Using 
McKenzie's  formulation,  the  range  of  correlations  obtainable 
with  the  EMA(l)  model  is  (— • 

An  alternative  procedure  for  producing  negative  correla¬ 
tions  was  introduced  by  Lawrance  and  Lewis  [Ref.  g] •  Their 
procedure  requires  two  series  of  innovative  factors  {En, 
n  =  0,1,...}  and  {E^,  n  =  0,1,...}  which  are  correlated  and 
may  be  antithetic. 

Antithetic  variables  are  generated  by  using  the  fact  that 
the  variables  E^  and  E|  defined  as: 


Ei  =  -ln(Ui), 

Ej  =  -InU-lL), 


(III.B.2) 


where  {U^,  i  =  1,2,...}  is  a  sequence  of  uniform  (0,1)  random 
variables,  are  both  marginally  Exponentially  distributed  and 
correlated.  P.A.P.  Moran  [3ef.  25]  has  determined  that  the 
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correlation  between  and  E|  is  approximately  -0.6449 

and  this  is  the  lowest  correlation  possible  between 

Exponential  random  variables.  E.  and  E?  have  a  deterministic 

i  1 

relationship  since 

-E. 

Ej  =  -In (1  -  e  x)  .  (III.B.3) 

Using  the  Lawrance  and  Lewis  extension  of  the  EMA(l) 
process ,  the  model  becomes 


X„  =  SE  +  I  E '  .  ,  (III.B.4) 

n  n  n  n-1 

where  {E^,  n  =  1,2,...}  is  an  iid  sequence  of  Exponential 

random  variables,  {E',  n  =  0,1,...}  is  a  sequence  of  Exponen- 

n 

tial  random  variables  which  are  correlated  with  the  respective 
variables  in  the  {En}  sequence,  {ln,  n  =  1,2,...}  is  a  sequence 
of  independent  binary  random  variables  with  P(ln=0)  =  1  -  P(In=l) 
=  3,  0  <_  S  <_  1,  and  {ln>,  {E^}  are  independent  of  each  other 
and  all  previous  Xn  values. 

The  correlation  of  the  X's  can  then  be  computed  as  follows. 

Let  E (X)  =  y  and  recall  that  since  {X_}  and  {E  }  have  the 

n  n 

same  marginal  exponential  distributions,  VAR(Xn)  =  VAR(En) . 


x„+ixn  - 


=  3  E  ±1£  +31  ..E.E'+BI  E  E  ’  7 +1  ^.1  E'E'  , 

n+1  n  n +1  n  n  n  n+1  n-1  n+1  n  n  n-1 
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Thus,  using  the  independence  of  {En},  {ln),  and  the  iid  nature 

of  (E  }  and  { I  } 
n  n 

E(Xn+lXn}  =  p2y2+3(l-3) [COV(En,E^)+y2]+B(l-6)u2+(l-3)2y2 

=  y2  +  3 (l-65COV(En,E^) 


Therefore , 


C0V(Xn+l'Xn)  =  3 (l-S)COV(En,E^) 

and  using  VAR(X)  =  VAR(E) 

CORR(Xn+1,Xn)  =  3 (l-B)CORR(En,E^)  (III.B.5) 

Using  Moran's  result  [Ref.  25]  the  correlation  of  antithetic 

Exponentials  is  known  to  be  approximately  -0.6449.  Therefore, 

the  greatest  negative  value  that  can  be  achieved  for  C0RR(Xn+^ ,Xn) 

is  approximately  -0.1612  when  3  =  0.5.  Since  no  restriction 

was  placed  on  C0RR(E  ,E'),  the  seauences  {E  }  and  {E1}  need 

n  n  n  n 

not  be  antithetic,  but  can  have  any  correlation  that  is  possi¬ 
ble  for  two  Exponential  sequences  with  the  same  marginal  dis¬ 
tribution.  By  specifying  that  E^  =  E^,  the  original  EMA(l) 
is  achieved  and  the  correlation  for  the  X's  is  3(1-3)  as  in 
Lawrance  and  Lewis  [Ref.  5],  By  allowing  the  correlation  be¬ 
tween  the  { }  and  {E^}  sequences  to  vary  from  -0.6449  to 
1.0,  the  correlation  of  the  X's  will  vary  from  a  minimum  of 
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-0.1612  to  a  maximum  of  0.25  (also  depending  on  the  value  of 
6)  as  can  be  seen  from  III.B.5.  The  Lawrance  and  Lewis  exten¬ 
sion  of  the  EMA(l)  model  gives  greater  possible  variation  in 
the  correlation  than  that  of  McKenzie,  but  requires  two  se¬ 
quences  of  Exponential  random  variables  to  achieve  this  range. 

Although  it  is  clear  that  with  E*  =  E  ,  a  C0RR(E  ,E')  =1 

n  n  n  n 

and  when  (E  }  and  { E 1 }  are  antithetic  CORR(E  ,E')  =  -0.6449, 
n  n  n  n 

it  may  not  be  obvious  how  to  generate  {E^}  sequences  with 
correlations  between  these  two  extreme  values.  A  simple  bi¬ 
variate  exponential  sequence  with  any  desired  correlation  in 
the  permissible  range  can  be  generated  using  the  relationship 

/  E^  with  probability  p, 

E|  =  !  (III.B.6) 

'  E^  with  probability  1-p, 

where  E^  is  the  antithetic  of  E^.  Then  the  extended  EMA(l) 
model  has  two  parameters,  3  and  p,  and  the  range  of  correla¬ 
tions  for  the  X's  is  -0.1612  to  0.25,  as  above.  The  bivariate 
density  for  the  pair  {E|,E^}  is  not  smooth.  Other  bivariate 
densities  such  as  those  in  Gaver  [Ref.  26]  and  Lawrance  and 
Lewis  [Ref.  27]  can  also  be  used. 

The  above  ideas  on  extending  the  correlation  structure  of 
the  moving  average  models  to  encompass  negative  correlation 
can  be  applied  to  all  of  the  new  models  given  below.  Details 
are  not  given. 
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C.  THE  EXTENDED  EMA(l)  MODEL,  NEMA(l) 

1 .  Introduction 

The  original  Exponential  moving  average  process  is 
discussed  by  Lawrance  and  Lewis  [Ref.  5] .  This  paper 
considers  the  first  order  process  defined  by: 

I  $En  with  probability  6, 

Xn  =  j  (III.C.l) 

(  3En  +  En+^  with  probability  (1-6), 

where  {En,  n  =  0,±1,±2,...}  is  an  iid  Exponential  sequence  and 

0  <_  6  <_  1 .  This  random  linear  combination  of  Exponential 

variates  is  called  the  EMA(l)  model  for  Exponential  moving 

average  of  order  one.  Since  X  is  a  function  of  E  and  E  , 

n  n  n+1 

this  version  is  called  the  forward  EMA(l) .  The  backward 
version  of  EMA(l)  is  obtained  when  En+1  is  replaced  by  En_^ 
in  III.C.l. 

The  fact  that  EMA(l)  is  a  single  parameter  model  sug¬ 
gests  that  this  model  may  not  be  sufficiently  flexible  to 
address  all  processes.  Investigation  of  an  alternate,  two 
parameter  model  may  indicate  that  a  two  parameter  model  is 
sufficiently  more  flexible  to  justify  its  increased  complexity. 

The  extended,  two  parameter  model  is  based  on  the  new 
Exponential  autoregressive  process  of  order  one  (NEAR(l) ) . 

The  NEAR ( 1 )  model  propounded  by  Lawrance  and  Lewis  [Ref.  8] 
is  defined  as 
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/  £n  +  SXn_^  with  probability  a 
Xn  =  j  (III. C. 1.2) 

'  En  with  probability  (1-a) 

where  (Xn,  n  =  1,2,...}  is  a  second-order  stationary  sequence 
of  Exponential  random  variables  with  parameter  X,  {En>  is  an 
iid  sequence  of  innovative  factors,  0  <  S  <  1,  0  <  a  <_  1,  and 
a+6  <2.  By  letting  <f>v(s)  and  <j>  (s)  denote  the  Laplace- 

«■  E 

Stieltjis  transform  of  X  and  E  respectively,  Lawrance  and  Lewis 
determined  that  <f>£(s)  =  (1-a'}  Bs  ‘  Usin9  a  partial 

fraction  solution  technique  to  invert  this  transform  produced 


i  E 

_  I  n 


with  probability 


(1-a 


(III. C. 1.3) 


(l-a)3En  with  probability 


where  {En,  n  =  0,1,2,...}  is  an  iid  sequence  of  Exponentially 
distributed  random  variables  which  has  the  same  parameter  as 
the  {Xn>  sequence. 

By  noticing  that  the  autoregressive  model  given  in 
III. C. 1.2  using  III. C. 1,3  is  a  random  sum  of  two  iid  Exponen¬ 
tially  distributed  random  variables,  the  NEMA(l)  model  is 
produced  by  substituting  En_^  for  Xn_-^  in  the  NEAR(l)  model. 
This  procedure  is  identical  to  that  used  to  produce  the  EMA(l) 
model  from  the  EAR ( 1 )  model  and  yields 


j  En  +  3En-l  with  probability  ot, 

X„  =  (III. C. 1.4) 

n  I 

’  En  with  probability  1-a. 
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This  model  can  be  written  in  a  more  compact  form  as 

Xn  "  KnEn  +  Wl'  (III. C. 1.5) 

where  {X  ,  n  =  1,2,...}  is  a  second  order  stationary  sequence 

of  marginally  Exponentially  distributed  random  variables; 

(En,  n  =  0,1,...}  is  an  iid  sequence  of  Exponential  random 

variables  with  the  same  parameter  as  the  {Xn}  sequence; 

{ln,  n  =  1,2,...}  is  an  iid  sequence  of  random  variables  with 

P  (I  =0)  =  1-P  (1=0)  =  a;  (K  ,  n  =  1,2, ...  }  is  an  iid  se¬ 
ll  n  n 

quence  of  random  variables  with  P(K  =  1)  =  1-P(K  =  (l-a)B)  = 

n  n 

8 

l-'(l-cir5';  ^n^'  ^Kn^'  and  ^En^  are  mutually  independent; 

0  <  d  <  1;  0  <_  0  <_  1;  and  a+0  <  2. 

The  NEMA(l)  model  contains  both  the  forward  and  back¬ 
ward  versions  of  EMA(l)  as  special  cases.  When  a  -  1; 

P(I  =0)  =  1,  (l-a)0  =  0,  and  P(Kn=0)  =  0.  Hence,  the 
NEMA(l)  model  becomes 

,  3En-1  with  probability  0, 

Xn  =  (III. C. 1.6) 

(  3En_^  +  En  with  probability  (1-0) . 

This  is  a  form  of  the  forward  EMA(l). 

When  0  =  1;  P(In=l)  =  a,  (l-a)6  =  (1-a)  ,  and 
P(Kn=  (1-a))  =  -g-  =1.  In  this  case,  the  NEMA(l)  model 

becomes 
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with  probability  (1-a)  , 


(l-a)E 


X 


( 

^  (l-a)En  +  with  probability  a. 


(III. C.  1.7) 


This  is  a  form  of  the  backward  EMA(l)  with  3  =  1-a.  There¬ 
fore,  the  NEMA(l)  contains  the  special  cases  of  EMA(l)  when 
a  and  8  assume  specific  values. 

One  should  also  note  in  passing  that  the  (Xn>  sequence 
becomes  an  iid  sequence  if  a  =  0  or  8  =  0. 

2 .  Correlation  Structure 

The  following  relationships  will  prove  of  value  in 
succeeding  calculations 


P(In=  8)  =  1  -  P(In=  0)  =  a. 
E(In)  «  a8  +  (l-a)*0  =  a8  • 

P(Kn=l)  =  l-P(Kn=  (l-a)B) 

=  1-3 

1- ( 1-a) 8  ' 

E‘V  “  + 

l-8+aS2-a232 
1-S+aS  " 


(III. C. 2.1) 

(III. C. 2. 2) 

(III. C. 2. 3) 


l-3+a3-a3+a82-a23'i 

I-8+a3 


l-8+a3  _  a3 (l-8+a8) 
1-8+aB  1-S+aB 
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E(K  )  =  1-aB  (III. C. 2. 4) 

Xn  "  KnEn  +  Wl  (III.C.2.5) 

E(Xn>  =  E(KnEn  +  Vn-l1  (III.C.2.6) 

=  E(KnEn)  +  Ed^^) 

=  E  (K  )  E  (E  )  +  E  (I  )  E  (E  ,) 
n  n  n  n- 1 

=  (l-a£) E (E  )  +  a3E (E  .) 

n  n-l 

E (X)  =  E(E) 

Since  {Xn;  and  {En:  are  both  Exponential,  E(X)  =  E (E)  implies 

VAR  ( X)  *  VAR(E) .  Since  both  (El  and  {X  }  have  the  same 

n  n 

Exponential  parameter,  without  loss  of  generality  this  param¬ 
eter  will  be  considered  to  be  1.  This,  of  course,  requires 
E(X)  =1  and  VAR{X)  =  1. 

The  possible  range  of  correlations  for  the  NEMA(l) 
model  can  be  determined  by  a  simple  calculation.  We  have 
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Using  the  independence  of  {Kn>,  {ln>  and  (En>  and  the  iid 
nature  of  these  three  sequences ,  we  have 

E(XnXn+i)  =  (l-a0)2[E(En) )2  +  (l-aB)aB [E(En) ]2 

+  (l-aB)aB [2{E(En) }2]  +  (ct3 ) 2 [E (En) ] 2 
=  1  +  (1-ad) ad 


Therefore, 

COV(Xn,Xn+1)  =  (l-ctB)aB 

and 

CORR (X  , X  , )  =  (l-aS)aB.  (III. C. 2. 7) 

n  n+1 

Therefore,  the  original  NEMA(l)  model  has  the  same  range 
of  possible  correlations  as  the  EMA(l),  namely  that  the  corre¬ 
lations  must  lie  in  the  interval  [0,^].  One  should  note  that  it 
is  not  possible  to  distinguish  the  parameters  a  and  3  from  the 
correlation,  or  even  whether  the  product,  ad,  has  a  given  value 
or  one  minus  that  value.  This  is  similar  to  the  normal  moving 


average  model  of  order  one  where  the  cases  <p  and  l/<p  are 

indistinguishable.  In  the  normal  case,  the  range  of  <J>  is 

limited  to  the  interval  [0,1]  on  the  basis  of  invertibility 

[Ref.  28] .  It  would  seem  simple  and  convenient  here  to  limit 

a8  to  the  interval  [0,i].  However,  non-normal  processes  are 

not  completely  determined  by  their  correlation  structure.  In 

fact,  Jacobs  and  Lewis  [Ref.  6  ]  showed  that  in  the  EMA(l) 

case,  the  values  8  and  (1-8)  can  be  distinguished  using  direc- 

2  2 

tional  moments,  E(xnxn+]_)  and  E^XnXn+l^‘  Hence»  such  a  re¬ 
striction  on  the  value  of  a6  is  inappropriate. 

One  should  also  note  that,  since  the  correlation 
between  Xn  and  Xn+K  is  zero  for  all  K  with  absolute  value 
greater  than  one,  the  first-order  correlation  completely 
determines  the  correlation  structure  of  the  model. 

The  restriction  on  the  range  of  attainable  correlation 
is  disappointing  but  not  surprising  since  it  can  be  proven  that 
any  Exponential  moving  average  process  of  order  one  generated 
as  a  linear  combination  of  independent  Exponentials  must  have 
a  correlation  that  lies  in  the  range  [0,^-].  The  proof  of 
this  contention  follows  the  previous  calculation  closely. 

THEOREM: 

Assume  {E  ,  n  =  0,1,2,...}  is  a  sequence  of  iid  Exponen¬ 
tial  variables  with  unit  mean, {An,  n  =  1,2,...}  and  (Bn, 
n  =  1,2,...}  are  sequences  of  iid  random  variables,  and  (An), 
{Bn},  (En }  are  all  mutually  independent.  Moreover,  assume  the 
sequence  {Xn,  n  =  1,2,...}  defined  by 


I 


Xn  =  A  E  +  B  E  . 
n  n  n  n  n- 1 


(III. C. 2. 8) 


is  a  unit  ir.ean  Exponential  sequence.  Now 


ElXn>  '  E<AnEn  +  Br,Vl> 


=  E  (A  )  E  (E  )  +  E (B  )  E (E  ,) 
n  n  n  n—  i 


1  =  E(An)  +  E(Bn) 


(III. C. 2. 9) 


In  addition,  since  XR  >  0  for  all  n,  both  An  and  Bn  must  be 

non- negative  for  all  n.  Hence  E(A  )  >0  and  E (B  )  >  0 .  It 

n  —  n  — 

also  follows  that  1  >  E(A)  and  1  >  E(B).  Now 


Therefore 


E(VW 


EfA^.A  E  ^tE  +  A  , . B  E  , E  .  +  A  B  ,  . E 
n+1  n  n+1  n  n+1  n  n+1  n-1  n  n+1  n 


B  .BEE  ,  ) 
n+1  n  n  n-1' 


[E  ( A)  ]  2  +  E  (A)  E  (B )  +  2E  (A)  E  (B)  +  [E  (B)  ]  2 


=  [E  (A)  +  E  (B)  P  +  E  ( A)  E  (B ) 


Since  E (A)  +  E(B)  =  1  from  III. C. 2. 9,  E (A)  +  E(B)  =  1, 


E(Xn  W  *  1  +  E  ( A)  E  (B) 
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Therefore,  since  E(XR)  is  one,  by  assumption, 

COV(Xn,Xn+1)  =  E (A) E (B) 

=  E  (A)  [1  -  E  (A)  J 

and 

CORR(Xn,Xn+1)  =  E ( A)  [ 1  -  E (A) ]  (III. C. 2. 10) 

Since  0  <_  E(A)  ^  1,  the  correlation  must  lie  in  the  interval 
[0,£]  .  Q.E.D. 

The  possible  range  of  correlations  can  be  extended  by 
reformating  the  model.  We  choose  to  do  this  first  by  using 
the  method  devised  by  Lawrance  and  Lewis  [Ref.  8]  and  given 
in  equation  III.B.4.  With  variables  and  sequences  defined  as 
in  equation  III. C. 1.5  let  {En>  and  {E^}  be  correlated  sequences 
and  redefine  the  NEMA(l)  as 

Xn  "  KnEn  +  VA-l*  .  (III. C. 2. 11) 

Then  it  follows  that 
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Thus 


E(XnVl> 


=  (1-ctB) 


+  (l-aB)aB  +  (l-aB)aB [COV [E„ ,E ’ ] +1] 

n  n 


2a  2 
a  3 


=  1  +  (l-a3)aBC0V(En/E^) 


and 


COV(Xn,Xn+1)  =  (l-a8)a3COV(En,E^) 


Fina  ly 


CORR(Xn,Xn+1)  =  (l-aB)otBCORR(En,E^)  .  (IIIC.2.12) 

As  described  above,  Moran  [Ref.  25]  has  determined  that 

the  range  of  possible  correlations  for  two  Exponentials  is 

(-0.6449,1.0).  Thus  when  aB  =  0.5,  the  possible  correlations 

for  X  and  X  fall  in  the  interval  (-0.1612,0.25).  This 
n  n+ i. 

procedure  extends  the  range  of  possible  correlations  at  the 
cost  of  generating  the  additional  {E^}  sequence. 

McKenzie  [Ref.  21]  has  suggested  that  the  range  of 
correlations  could  be  extended  by  requiring  that  the  {I  } 
sequence  be  correlated.  Using  this  scheme,  he  was  able  to 
show  that  the  correlations  for  the  (Xn>  sequence  lies  in  the 
interval  (-  •  Because  of  the  requirement  of  the  moving 

average  process  of  order  one  to  have  zero  correlation  for  lags 
greater  than  one,  the  correlation  of  the  { In }  sequence  also 
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had  to  be  of  MA(1)  type.  It  is  this  restriction  that  pro¬ 
duces  such  a  narrow  range  of  correlations.  A  logical  and 
obvious  extension  to  McKenzie's  work  is  to  require  that  both 

the  {K  }  and  {I  }  sequences  in  the  NEMA(l)  model  have  a  MA(1) 
n  n 

correlation  structure.  This  can  be  combined  with  the  corre¬ 
lated  {E  },  { E ' }  scheme  of  Lawrance  and  Lewis.  If  this  com- 
n  n 

bined  case  is  carried  out,  the  NEMA(l)  model  is  as  follows 

Xn  -  Vn  +  VA-l'  (III. C. 2. 13) 

where  {Kn,  n  =  1,2,...}  is  a  sequence  of  random  variables 
with  an  MA(1)  correlation  structure  with  P(Kn=l)  = 

l-P(Kn=  ( 1-a)  3 )  =  XMT^aT8;  {ln'  n  =  is  a  sequence 

of  random  variables  with  an  MA(1)  correlation  structure  with 

P(I  =  B)  =  1-P (I  =  0 )  =  a;  {E  }  and  { E *  >  are  correlated  se- 
n  n  n  n 

quences  with  marginal  Exponential  distributions  with  unit 
means;  and  {Kn},  {1^},  and  { Er }  are  mutually  independent. 

Now 
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E(XnXn+l}  =  E(Kn+1Kn)+(l-aS)aS+(l-a6)a6E(EnE^)+E(In+1In) 

=  COV(Kn+1,Kn)+(l-aB) 2+(l-a6)ag 

+  (l-aB)a6-  [COV (En ,E^) +1] +COV (I  ,  In>  +  (aB) 2 

=  l+COV(Kn+1,Kn)+COV(In+1,In)+(l-aB)aBCOV(En,E^) 

Therefore 

COV(Xn,Xn+1)  =  COV(Kn+1/Kn)+COV(In+1,In)+(l-aB)oBCOV(EnfE^) 

and 

CORR(Xn,Xn+1)  =  COV(Kn+1,Kn)+COV(In+1,In)  (III. C. 2. 14) 

+(l-aB)aBCOV(En,E^) 

Although  this  scheme  obviously  extends  the  range  of  possible 
correlations,  it  does  so  at  the  expense  of  considerable  com¬ 
plexity.  Considering  the  limited  range  of  correlations 
possible  by  imposing  a  correlation  on  the  {I  }  and  {Kn> 
sequences,  the  additional  complexity  may  not  be  warranted. 

If  in  spite  of  the  complexities  involved,  one  decides  to  in¬ 
duce  correlations  in  the  coefficient  sequences,  the  NEMA(l) 
because  it  has  two  such  sequences  will  yield  a  slightly  larger 
range  than  the  EMA(l)  model. 
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3. 


?(Vl>Xn! 

One  of  the  possible  advantages  of  a  two  parameter 
model  is  the  capacity  for  modifying  P(Xn+^>Xn)  and,  conse¬ 
quently,  the  sample  path  behavior  of  the  process  while  main¬ 
taining  a  constant  correlation.  Since  the  correlation  is  a 
function  of  a6 ,  one  can  vary  the  values  of  a  and  3  while 
keeping  the  product  constant.  The  P (Xn+^ > XR)  can  be  calcu¬ 
lated  by  addressing  each  of  the  sixteen  possible  combinations 
of  K  and  I  values  for  Xn+^  and  Xn,  computing  the  probability 
for  each  combination ,  and  weighting  the  probability  associated 
with  a  given  combination  by  the  probability  that  the  given 
combination  occurs.  A  sample  calculation  is  provided  and 
complete  results  presented  in  Table  III. C. 3.1. 

Example:  We  have 

X  =  K  E  +  I  E  ,  , 
n  n  n  n  n- 1 ' 


Xn+1  * 

Kn+lEn+l  + 

1  .  i  E  , 
n+1  n 

P(Kn=D 

-  X-p(Kn 

=  ( 1-a) 6 )  = 

1-6 

1- (1-a) 6  ' 

P(In“  3) 

“  X-P(1n 

*  0)  =  a. 

Let 


1-6 

T-TL-^6 


=  '5  and  consider  the  case  where 


Xn  "  3'  Kn  =  lf  In+1  3'  Kn+1  L* 
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Since  the  {I  }  and  {K  1  sequences  are  both  iid  and 
n  n 

independent  of  each  other,  the  probability  of  this  combination 
of  parameter  values  is  simply  the  product  of  the  individual 
probabilities  of  occurrence.  Hence  the  probability  of 
occurrence  is  a  5  .  Then  in  this  case 


*<*„+! 'V  *  p<En+l+6En>En+SEn-l> 


-  p<En+l>  a-i)EntaVl' 


(III. C. 3.1) 


Now  En+-^  is  independent  of  Y  =  (1-3)  En+^n-i  •  Therefore,  the 

calculation  required  by  equation  III. C. 3.1  is  straightforward 

once  the  density  of  Y  is  obtained.  We  have,  with  fE  (x) 

n-1 

the  p.d.f.  of  En-1  (i.e.  e  X) 


y/3 


P(  [1-B]E+3E  ,  <_  y)  =  /  P(  [1-3]E  +3x<_yiEn-1  *  x)f£  (x)dx 

n  n  i  g  n-1 


y/3 


=  /  P([1-3]E  <y-3x;En_i  =  x)fE  (x)dx 

n  n-1 


y/3 


p  ( e  <  Xfl^iEn  .  =  x)  fF 

;0  n-  1-b  1  n-1  En_i 


(x)  dx 


,,  y-3x 

y/3  -x 

/  ( 1  -  e  1  )  e  dx 

'  0 

y/3  . 


y/3 

e  xdx  -  /  e  “e  dx 

0  0 


y  ,,  r  1-23 1 


=  1  -  e 


-  e-y/3  - 


j  e 
0 


dx 
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P( [l-8]En+BEn_1 <  y) 


=  1-e-^8- 

=  l-e^8- 

=  l-e-^8- 

P  (Y  <  y)  =  - 


,  rY/e  rl-2B  i  _X[T^'] 

*  lt^b]  >q  hrr16 


;A,  i- 


3  r  1~  8  i  i  1 
[ i- 2 6  (1~e 


_Zr1-2B1 

6[T^T] 


,^,e-A+ 


4^[l-e~y/S]  + 


1^21 


dx 


Therefore,  the  density  function  of  Y,  f^Cy),  is 
-  ( 3 )  (|pe  ^ X' ~~3  e  ^ 1  ^  »  for  8  ^  Gaver  and 

Lewis  [Ref.  2]  gave  the  necessary  and  sufficient  conditions 

“ A . X  - A_X 

for  a  mixed  exponential  of  the  form  TT^A^e  +’T2^2e  ' 

>  1,  tt^  +  rr ^  =  1,  and  <  A2  to  be  a  proper  density 
function.  The  condition  is  that 


*1  - 


(III. C. 3. 2) 


In  this  situation  we  address  two  separate  cases.  The  first 
case  is  when  0  <  3  <  In  this  case  r—- ^  >  1  and  the 

requirement  is 


-1 


[1  - 


-1 


,  1-3 

-  fib 
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p(xn+l>Xn)  ~  '(rTF)(BTT)  +  lI^J)  i^) 

"Vl’V  -  (S-M1U6)  (III. C. 3. 3) 

Table  III. C. 3.1  presents  the  results  of  the  calculations  for 

all  of  the  sixteen  combinations  of  parameter  values  for  Xr+1 

and  X  . 
n 

When  the  P(X  >  X  )  was  calculated  for  various  values 
n+i  n 

of  a  and  0,  it  was  found  that  the  values  for  this  probability 
varied  from  0.44  to  0.54.  Table  III  .C  .3 . 2  contains  the  results 
of  these  calculations  for  four  hundred  forty-one  combinations 
of  parameter  values.  Although  the  variation  in  probability  is 
not  large,  it  does  represent  an  increase  over  the  forward 
EMA(l)  model.  In  particular,  the  forward  EMA(l)  model  can  not 
produce  a  probability  greater  than  0.50.  Consequently,  the 
NEMA(l)  model  not  only  has  a  greater  range  of  possible  proba¬ 
bilities,  but  also  can  produce  probabilities  greater  than  0.50. 
The  implications  of  this  greater  range  is  that  the  NEMA(l)  model 
can  address  data  sample  paths  that  have  a  slight  tendency  for 
either  runs  of  increasing  or  decreasing  values,  while  the 
EMA(l)  can  only  address  sample  paths  that  tend  to  produce  runs 
of  decreasing  values. 

Examples  of  scatter  plots  and  sample  paths  for  three 
sets  of  parameter  values  and  positive  correlations  are  given 
in  Figures  III .C . 3 . 1-III.C . 3 . 6 .  Because  of  the  relatively  low 
correlations  possible  and  because  of  the  limited  range  of 
values  for  P(Xn+1>Xn),  differences  among  the  figures  are  not 
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sharply  delineated.  However,  differences  can  be  detected, 

particularly  in  the  sample  paths.  In  Figure  III. C. 3.1  with 

an  a  value  of  0.95  and  S  value  of  0.50  has  a  P(X  >  X  )  = 

n+1  n 

0.44,  the  lowest  value  for  this  probability.  Since  this  pro¬ 
duces  a  slight  tendency  for  runs  of  decreasing  value,  the 
number  of  extreme  values  (i.e.  greater  than  3.0)  is  two. 

In  Figures  III. C. 3.2  and  IIIC.3.3  the  P (X  ^ >  Xn)  is  0.50  and 
0.54,  respectively,  with  a  corresponding  increase  in 
the  number  of  large  values.  This  trend  is  more 
difficult  to  detect  in  the  corresponding  scatter  plots. 
Figures  III .C . 3 . 7-III .C  .  3 . 12  provide  sample  paths  and  scatter 
plots  for  the  same  a  and  3  values  as  previously  displayed, 
but  with  antithetic  innovative  sequences  (see  III.B)  and 
consequent  negative  correlations.  Although  the  negative 
correlation  is  evident,  trends  in  these  figures  are  difficult 
to  detect.  The  extremes  of  sample  path  variability  produced 
by  the  NEAR ( 1 )  process  [Ref.  8]  are  not  reproducible  with  the 
NEMA(l)  process.  This  may  be  attributable  to  the  restricted 
range  of  possible  correlations. 

4 .  Laplace  Transform  of  Sums 

One  aspect  of  the  EMA(l)  model  is  its  analytical 
tractability .  This  is  evidenced  by  the  ability  to  derive  the 
Laplace  transform  of  sums  by  a  recursive  relationship  given 
by  Lawrance  and  Lewis  [Ref.  5] .  This  tractability  carries 
over  to  the  NEMA(l)  process.  The  Laplace  transform  is  useful 
in  obtaining  quantities  which  are  of  use  in  analyzing  point 
processes,  namely  the  intensity  function  and  the  spectrum  of 
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FIGURE  III. 
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SRMPLE  RHO:  0.26 
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X  (N) 


FIGURE  III. C. 3. 4 
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FIGURE  III. 


FIGURE  III. 


FIGURE  III. C. 3. 11 
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FIGURE  II. C. 3. 12 
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counts.  These  quantities  are  derived  for  the  NEMA(I)  process 

in  subsequent  sections  from  the  results  obtained  here. 

1  —  s 

With  X  defined  as  in  III. C.  1.5,  let  ■. — n — mr  =  3 

n  r  l-(l-a)3  _sT 

and  (l-a)8  =  y.  Further,  let  £  X.  =  T  and  (s)  =  E(e  r)  . 

i=l  1  r  r 

Then  we  have 


T  =  X.  +  X-  +  .  .  .  +  X 
r  1  2  r 


(III. C. 4 .1) 


K1E1  +  IjEq  +  K2E2  +  I2E1  +  ...  +  KrEr  +  1^^ 


KrEr  +  (Ir  +  Kr-l)Er-l  +  •••  +  d2  +  Ki)El  +  IiE0 


Then  letting  =  Ij  +  ^  +  Kj/  i  =  1/2,..., r-1  and  using  the 
mutual  independence  of  the  iid  sequences  { } ,  {ln},  { En } 


—  sT 

4>t  (s)  *  E(e  r) 


(III. C. 4. 2) 


-S[K  E  +L  .E  ,  +  .  .  .+L.E.+I.E.] 

E  (e  r  r  r~L  1  ) 


-sK  E  ~sEr-lEr-l  -sL.E^  -sI^Eq 

=  E  (e  r  )E(e  r  1  i)...E(e  1  )E(e  1  ) 


Now  let  VKis)  =  E(e 
Then 


-SK.E. 
3  3 


-si -E . 
3  3 


) ,  (s)  =  E(e  J  J) ,  V. (s)  =  E(e 


-sL.E . 
3  3 


$T  (s)  =  ^(s)  f x  ^L(s)  ^ 


r-1 


(III. C. 4. 3) 


To  evaluate  these  quantities  note  that 
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I 


I 


K 


j  1  with  probability  5 , 

^  Y  with  probability  1-6 . 


Then 


¥k(s)  =  E(e 


-sVi, 


-SE.  -SYE. 

=  6E(e  :)  +  (1-6 ) E (e  J ) 


4*  ( s )  =  6$„(s)  +  (1-6)  <p  (ys)  , 


K 


where  <frE(s)  =  j+g 


.  So 


Vs) 


<S  +  (1-6) 
1+S  1+Y  S 


Also 


(III. c. 4. 4) 


(III. C. 4. 5) 


|  3  with  probability  a, 

1  =  / 

'  o  with  probability  1-a, 

so  that 

(III. C. 4. 6 

(III.  C. 4. 7) 


-sI.E.)  -s8E. 

4/^(3)  =  E(e  ^  ^  =  aE(e  J)  +  (1-a) 


'tI  (s) 


•s — 2L —  +  (l-a) 
l  +  3s 
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To  evaluate  'V-  (s)  note  that 

■U 


/  3+1  with  probability  a6 , 

1  3+y  with  probability  ad-6), 

L  *  \ 

I  1  with  probability  (1-a) 6 , 

Y  with  probability  (1-a)  (1-6)  . 

Therefore 

-SL.E.) 

i'L(s)  =  E(e  33 

-s [6  +  1] E  .  -s [6+Y 1 E , 

=  a6E(e  3 )  +  a(l-6)E(e  3  ) 

-SE.  -syE. 

+  (1-a) 6E (e  3)  +  (1-a) (l-o)E(e  3 ) 

«L(s)  =  a5<DE([B+l]s)  +a(l-5)0E([3+Y]s)  +  (l-a)5DE(s)  (III. C. 4. 8) 

+  (1-a)  (1-6)  4>e(ys)  . 

Using  the  results  of  III. C. 4. 3,  III. C. 4. 4,  III. C. 4. 7,  and 
III. C. 4. 8 

6t  (s)  =  [6$e(s)  +  (1-6)4>e(ys)  ]*  [a<DE(3s)  +  (l-a)  ]  (III. C. 4. 9) 

1  r 

X [a6$E( [8+l)s)+a(l-6)$E( [0+Y]s)+(l-a)6*E(s) 

+  (1-a) (1-6)De(ys) ]r_1 
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This  extends  the  result  (3.7)  in  Lawrance  and  Lewis  [Ref.  5] 
to  the  NEMA(l)  process. 

5 .  Laplace  Transform  of  the  Distribution 
of  Counts 

The  Laplace  transform  of  the  sum  is  useful  in  deriving 
the  distribution  of  the  synchronous  counting  process  of  the 
number  of  events  that  occur  in  (0,t]  when  the  origin  is  estab¬ 
lished  at  the  occurrence  of  an  arbitrary  event.  The  number 
of  events  in  (0,t]  is  related  to  the  distribution  of  a  sum  by 
the  relationship 

<  r  iff  Tr  >  t,  r  =  1,2,...  (III. C. 5.1) 

f 

where  is  the  number  of  events  in  (0,t]  and  T^  is  the  sum 
of  the  first  interevent  times.  Thus 

P(N^=  r)  =  Fr(t)  -  Fr  +  1(t) 


where  Fr(t)  is  the  distribution  function  of  Tr>  The  proba¬ 
bility  generating  function  of  can  then  be  written  as 


=  l  ZrP(N^  =  r) 
r=0  c 


00 

=  l  Zr[F  (t)  -  Fr+1(t)] 
r=0 

00  l 

=  1  +  (Z-l)  l  Zr-1F  (t) . 

r=l  r 
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*  ★ 

Let  4'f(z;s)  be  the  Laplace  transform  of  4'f(z;t),  and  fr(t) 

the  Laplace  transform  of  fr(t),  the  p.d.f.  of  Tr>  Then 

•k  00 

^(z;s)  =  4  -  ■1~~)  l  zr-1f*(t)  (Hi. C. 5. 2) 

r  s  s  r=l  r 

Using  the  Laplace  transform  of  the  sum  from  III. C. 4. 9 

00 

•  £  ( z ;  s )  =  ■1~2-)  I  Zr_1  [d4>  (s)  +  (1-6  )  4>  e(YS)] 

r=l 

X  [ at 4) E  (Ss)  +  (1-a)  ]  x  [a6<j>E  (  [3+1]  s)  +a  (1-6  )  <j>E  (  [3+y  ]  s) 

+  (l-a)5  4>E(s}  +  (l-a)  (1-6 )  <J>E  (ys)  ]  r_1 


•k 

V  f  ( z ;  S ) 


i“  -~^-lfi4iE(s)  +  (l-5)0E(Ys)  Jx  [aq>E(Ss) 


(III. C. 5. 3) 


+  (l-a)  ]  xl 


)l-z  Id6<pE  (  [0+1]  s+a  (1- 
+  (l-a)e\,  (s)  +  { 1-a)  (1- 


—5 )  4>e  (  t3  +y )  s) 


6 )  <J>E (y s)  ] 


where  i)£(s)  = 

If  m^(t)  is  the  intensity  function  of  the  point  proc- 
★ 

ess,  then  mf(t),  its  Laplace  transform,  can  be  obtained  by 
differentiating  III. C. 5. 3  with  respect  to  z,  evaluating  the 
derivative  at  z  =  1,  and  then  differentiating  with  respect  to 
s.  These  steps,  when  taken,  produce  a  series  of  tedious 
calculations  which  produce  no  analytical  insights.  The  result 
of  these  steps  is 
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This  result  can  be  verified  in  a  number  of  ways.  First,  when 

3=1,  the  process  is  the  EMA(l)  process  and,  hence,  the  formula 

must  reduce  to  (4.2)  given  in  Lawrance  and  Lewis  [Ref.  5]  with 

X  =  1.  Second,  when  a  =  0,  the  NEMA(l)  process  reduces  to  a 

Poisson  process  and  the  formula  under  this  condition  must 

reduce  to  the  Laplace  transform  of  the  constant  intensity 

function  of  a  Poisson  process  with  rate  1,  j.  Third,  with 

3  =  0  the  NEMA(l)  process  is  again  a  Poisson  process.  Finally, 

* 

using  one  of  the  Tauberian  Theorems,  lim  mf(t)  =  lim  smf(s)  =  1. 

t-^oo  r  s->-0  r 

We  take  these  cases  in  turn.  First,  when  a  =  1 

l+(l+48-2a8)s+(3S+5B2-2aB-5a62+a2B2)s2 

+  (2B2+2S3-3aS2-3a33+ct232+a283)s3 _ 

s+(l+46-3a3+a232) s2+ < 38+58 2-2a6-7a32+3a2B2) s3 

2  3  2  322234 

+(2S  +23  -3aS  -3a6  +a  3  +a  3  ) s 


* 

mf(s)  = 
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* 

mf  (s) 


(1+68) ( 1+ [1+8 ] s) 
s ( [1+8] [1+ (8+82) s] ) 


(i+Bs)  (i+  [i+e j s) 

«(wl»llts"riraTtsl' 


which  is  the  result  of  Lawrance  and  Lewis  [Ref.  5]  with  A  *  1. 
In  the  second  case  with  a  =  0 


* 

m_  (s) 


l+(l+48) s+(3B+532)s2+(2B2+283)s3 
s+  (1+4S) s2+( 36+56 2)s3+(2B2+2B3)s4 


1 

s' 


the  Laplace  transform  of  a  Poisson  process  with  rate  of  1. 
In  the  third  case  6  =  0 ,  so 


mf  (s) 


1  +  s 
s  +  s4 


again  the  Laplace  transform  of  a  Poisson  process  with  a  rate 
of  1 . 

In  the  final  case  apply  the  Tauberian  Theorem 


l+(l+48-2a8) s+(3B+532-2aB-5aB2+a232)s2 


lim  smf (s) 
s-0 


2  3  222233 

+  (23  +26  -3a8  +a  8  +a  6  ) s _ 

l+{l+46-3aS+a232) s+ ( 38+5B2-2aB-7aB2+3a232 ) s2 


2  3  222233 

+  (2S  +23  -3ctS  +ct  3  +a  3  )  sJ 


1 

T' 


as  required. 
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6 .  The  Spectrum  of  Counts 

For  the  statistical  analysis  of  series  of  events  the 

most  useful  quantity  associated  with  a  process  is  the 

(Bartlett)  spectrum  of  counts.  The  spectrum  of  counts, 

g+ (w) ,  is  the  Fourier  transform  of  the  covariance  density  of 

N^(t).  It  is  related  to  the  Laplace  transform  of  the  inten- 
★ 

sity  function,  m^(s),  by  the  relationship  derived  by  Cox  and 
Lewis  [Ref.  29] 

\  *  * 

g+  (oj )  =  —  (1  +  m*  [iu>]  +  m^  [— iui ] )  . 


We  now  derive  this  for  the  NEMA(l)  process  using 

★ 

In  that  expression  for  m^(s),  let 
a.^  =  1  +  4S-2aS, 

a2  =  33  +  582  -  2a8“  5a32  +a232, 

a3  =  282  +  28 3  -  3aS2  -  3ct83  +  a282  +  ct2B3, 

bx  =  1  +  48  -  3aS  +  a232, 

b2  =  33  +  582  -  2ct8  -  7a82  +  3a282, 
b3  =  282  +  283  -  3a82  -  3a63  +  ct2S2  +  ct233. 


III. C. 5. 4. 


(III. C. 6.1) 

(III. C. 6. 2) 

(III. C. 6. 3) 

(III. C. 6. 4) 

(III. C. 6. 5) 


(III. C. 6. 6) 


Then 
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Recall 


<J+(w) 


g+  (w ) 


j 


*  1+a, s+a_s  +a,s 

_  ,_v  _  1  2 _ 3 

nu(s)  -  2  3  2 

s+b^s  +b2s  +b3s 


that  A  =  1  so 


2  3 

^  1+a^  (iw )  +a2  (iw )  +a3(iw) 

—  —  [  1  +  - * - S - T 

T  iw+b^iw)  +b2  (iw)  +b3  (iw) 

2  3 

1+a^  (-iw )  +a2  (-iw )  +a3(-iw) 

+  2  3  J- 

-iw+b1(-iw)  +b2(-iw)  +b3(-iw) 


2  3  4 

[iw+b^(iw)  +b2(iw)  +b3(iw)  ] 

2  3  4 

j  ><-iw+b1  (-iw)  +b2(-iw)  +b3(-iw)  ] 

'r  '  [iw+b1(iw)2+b2(iw)3+b_  (iw)4] 

|x  [-iw+b1  (-iw)  2+b2  (-iw)  3+b3  (-iw)4] 


2  3 

[1+a-^  (iw) +a2  (iw)  +a3(iw)  ] 

+  x  f-iw+b^  (-iw)  2+b2  (-iw)  3+b3  (-iw)  4  ] 
[iw+b^  (iw)  2+b2  (iw)  3+b3  (iw)  4  ] 

x  [  —  iw+b ^  (  —  iw )  2  +b 2  (  —  iw )  3+b 3  (  —  iw)  4] 


2  3 

[1+a^ (-iw) +a2 (-iw)  +a3 (-iw)  ] 

2  3  4 

x(iw+b^(iw)  +b2(iw)  +b3(iw)  ] 

2  '  3  4 

(iw+b1(iw)  +b2(iw)  +b3 (iw)  ] 

x  [-iw+bi (-iw) 2+b2 (-iw) 3+b3 (-iw) 4 ] 
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Consider  the  first  term  with  numerator  and  denominator  the 
same.  Let  the  denominator  =  D. 


D  =  [ioj+b^  ( icj )  2+b2  ( icu )  2+b^  (ico)  4]  [-iw+b^  (-iw)  2+b2  (-iw)  ^-t-b^  (-iu)  4] 


2  . ,  3  ,  4  5 ,  3 ,  ,2  4  ...  5  ,,  6  .  4 

gj  -ib^w  -b2GJ  +ib2<jj  +ib^Gj  +b^w  -ib^b2w  -b^b^uj  -b2GJ 


526  7  5  6  728 

+ib^b2GJ  +b2GJ  -ib2b.3GL>  -ib^w  -b-^b^oj  +ib2b2w  'H-b^w 


353  5  5  7  75 

i(-b^o)  +b2GJ  +b^o)  -b^b2<jJ  +b^b2oj  -b2b2Go  +b2b3W  ) 


2  ,  4  ,  ,2  4  .  ,  6  ,  4.2  6  .  .  6^,2  8, 

+  (gj  -b2oJ  +b^oj  -b^b^w  -b2uj  +^2W  ^31*1  ) 


gj2  [l+[b2-2b2]u2+[b2-2b1b3]aJ4+b2aJ6) 


(III.C.6.8) 


Let 


X  =  l+(b2-2b2)u)2+(b2-2b1b3)u)4+b2G)6, 

where  b^,  b2,  and  b^  are  defined  by  III. C. 6. 4,  III. C. 6. 5, 

III. C. 6. 6,  respectively.  Then 

D  =  u)2X  (III. C. 6. 9) 


Consider  the  numerator  of  the  second  term  in  III. C. 6. 7  and 
call  it  N2.  Then 
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N2  =  [1+a^  (iaj) +a2  (iui)  (ioj)  [-iw+b^  (-iw)  ^  +  b2  (~iw)  (-iuj)  4] 


2342  3  4  53 

-lu-b^uj  +^2^  +b2^  +aiUJ  -ia^b^oj  -a^^oj  tia^^on  +ia2^ 

+a2b^a)^ -ia2b2W^-a2b a^m^+ia  ^^(^+3^2^“  ia^b  2^ 


N2  =  i  {-'^+ [a2~a^b^+b2]  (a^b2~a2b2+a2b^)  o)^-a2b2w'7)  (III.C.6.10) 


2  4  6 

+  (a^-b^)u  +  (a2b^-a^b2 -a2+b2 )  j  ^  > 


where  a^ ,  a2 ,  a b^,  b2,  and  b3  are  defined  in  III. C. 6.1 
through  III. C. 6. 6  respectively.  Consider  the  numerator  of 
the  third  term  in  III. C. 6. 7  and  call  it  N3.  Then 


N3  =  [1+a^  (-iuj) +a2  (~iw)  ^+a2  (~i^)  [iw+b^  (iui)  ^+b2  (iw)  ^+b2  (iui)  ^  1 


2342  3  4  53 

iuj-b^w  -ib2iJJ  +b2^  +a^w  +ia^b^w  -a-j^w  -ia.j^'-o  -ia2^ 

+a2b1Jjii+ia2b20i)^-a2b2^6-a2^4-ia2b^aj^+a2b2-J^+ia2b2^7 


N3  =  i  (-J- [a2-a1b1+b2]  w3- [a1b2-a2b2+a2b1]  uj5+a2b2W7)  (III.  C.  6. 11) 


2  4  6 

+  ( ai“bi )  w  +  (a2b^-a^b2- a-+b3)  J  + {a2b2-a2b2 ) w  , 


where  a^,  a2»  a,,  b^f  b2»  and  b2  are  defined  in  III. C. 6.1 
through  III. C. 6. 6,  respectively.  Note  from  III.  C.  6. 7  that 
all  terms  in  the  sum  have  the  same  denominator.  Use  III.C.6.10 
and  III. C. 6. 11  to  determine  the  numerator  of  the  second  and 
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third  terms  and  call  it  N.  Then 

3  5  7  3 

N  =  i  [-u)+  (a2_a^b^+b2)  w  + (a^b3~a2b2+a3b^ ) w  -a^b^ui  +  oo-  (a2~a^b^+b2)  w 

5  7  2  4 

-  (a1b3-a2b2+a3b1)  to  +a3b3w  ]+(a^-b^)w  +  (a2b1-a1b2-a3+b3)  oj 

+  (a3b2-a2b3)oj^+(a^-b^)w^+  (a2b^-a^b2-a3+t>3)  u)4+(a3b2-a2b3)cj6 

2  4  6 

=  2[(a^-b^)w  +  (a2b^-a^b2“a3+ti2 )  w  +  ^a3^2~a2b3 )  w  ^  • 

Let 

y  =  (a^-b^)  +  (a2b^-a1b2~a3+t>3  ^  (a3b2-a9i>3  ^ 01)6  ^  •  (HI.C.6.12) 

Then 

N  =  2oj2y  (III.  C.  6. 13) 


Using  III. C. 6. 7,  III. C. 6. 8,  III. C. 6. 13 


g+(w) 


.2  -2 

1  ,0J  X  ,  2^  v, 

If'  2  + 

a)  X  il)  X 


(III. C. 6. 14) 


where  x  and  y  are  defined  in  III.C.6.  9  and  III. C. 6. 12, 
respectively . 
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Figures  III .C. 6.1  through  III. C. 6. 3  show  the  results 
of  the  calculation  of  the  (Bartlett)  spectrum  of  counts.  In 
presenting  the  results  the  constant  j  in  III. C. 6. 14  was  ig¬ 
nored.  Figure  III. C. 6.1  shows  the  spectrum  of  counts  for  the 
same  a  and  3  values  that  were  used  for  the  sample  paths  and 
scatter  plots  of  Figures  III.  C.  3.1  through  IIIC.3.6.  This 
figure  also  shows  the  variation  in  the  spectrum  of  counts  as 

the  P(X„,T  >  X  )  varies  from  its  lowest  to  highest  values, 
n+l  n 

Figure  III. C. 3. 2  holds  the  P(X  .  >  X  )  constant  and  varies 
^  n+l  n 

the  correlation.  Since  the  spectrum  of  counts  for  a  Poisson 
process  is  a  constant  one  when  X  equals  1  and  the  constant 
—  is  ignored,  the  correlation  can  be  viewed  as  a  measure  of 
the  process '  departure  from  a  Poisson  process .  This 
divergence  as  a  function  of  the  correlation  shows  clearly  in 
this  figure.  Figure  III. C. 6. 3  holds  the  correlation  constant 
and  varies  the  P(Xn+1>  XR) .  The  slight  variation  in  the 
spectra  shows  that  while  the  spectrum  of  counts  does  vary 
with  the  P(Xn+^>  X  )  ,  the  correlation  plays  a  more  dominant 
role . 

The  analysis  from  the  Laplace  transform  of  sums  in 
III.C.4,  through  the  Laplace  transform  of  the  intensity  func¬ 
tion  in  III.C.5,  to  the  spectrum  of  counts  in  this  section 

» 

can  be  performed  using  the  correlated  {En),  (En)  se<3uences 
of  III.C.2  and  thus  for  negative  correlations.  Details  are 
not  given. 
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FIGURE  III 


i 


7 .  Joint  Laplace-Stieltjes  Transform  of  Xn  and  Xn|^ 


Because  the  NEMA(l)  process  is,  by  construction, 
only  one-dependent,  all  of  the  second-order  properties  of 
{Xn>  are  contained  in  adjacent  pairs  {Xn,Xn+^}.  In  previous 
sections  quantifiers  of  the  distribution  of  {Xn,Xn+^}  such 
as  p.  and  P(Xn  +  ^  >  xn^  have  been  derived.  Here  we  give  the 
Laplace-Stieltjes  transform  of  the  joint  distribution.  One 
could,  for  example,  study  the  effect  of  the  two  parameters 
from  this  result  by  deriving  directional  moments. 

The  joint  Laplace-Stieltjes  transform  of  Xn  and  Xn+1 
can  be  calculated  by  considering  each  of  the  sixteen  possible 
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1 


-sl[En+6En-l]"S2[En+l+6En] 

(s,  ,s.)  -  aa66E  (e 

Wl  1  2  ,  r 

+  <l-a)a«SE(e  ) 

„„  -A-!  W , 

+  a  (1-a)  <56E  (e  * 

*S-1  E  —  s ,  i 

+  (1-0) (l-a)S«E(e  1) 

,  -=1l'lV,En-l1'!2lVl,!V| 

+  aa6  (1-6  )  E  (e 

-sLlYVBEn-l]'S2En+l, 

+•  ( 1-a ) a6 ( 1-5 )  E  (e 

-«1TTEn-.2  tEn+l+6Enl 
+  a  (1-a)  6  (1-6 )  E  (e  ' 

“slYEn_s2En+l. 

+  (1-a)  (1-a)  6  (1-6  )  E  (e  ) 

-sl!V3En-l1"s2  lTEn-H+SEn) . 

+  oo(l-J)iE(e  1 

-sl[V6En-l1'Wl, 

+  (l-o)o(l-«)5E(«  1 

-»lE„-*2  "Vi*®.1, 

+  a(l-a) (l-5)5E(e  x  > 

_slEn_s2yEn+lv 
+  (1-a)  (1-a)  (1-6 ) 6E (e  > 

■sl[YEn+SEn-l1-S2[yEn+l+BEn]  . 

+  aa (1-6 ) (1-6 ) E (e 

.  -si[^En+BEn-l]'S2YEn+l) 

+  ( 1-a) a ( 1-5 ) (1-6 ) E (e 

■sl'Vs)1,Vlti£n1, 

+  a(l-a) (1-6) (l-5)E(e  > 

..  ,  ~ s  1Y En” S 2 Y En+ 1  > 

+  (1-a) (1-a) (1-6) (1-5) E (e  > 
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4>v  v  (s. /S-)  =  aa65ipE(s1+Bs2)  EgCSSj^)  eE(s2) 

Xn'  n+1  L 

+  (l-a)a654>E(s1)  <f)E(Ss1)  $E(s2) 

+  a  (1-a)  56  4>E  (s^+0s2)  4>E  (s2) 

+  (1-a)  (l-a)564»E(s1)«E(s2) 

+  aao  (1-5)  4>E  (YS2+8s2)  0E  (BSj^)  0E  (s2) 

+  (1-a)  a6  (1-6)  ^(ys^  4>e(6s1)  4>e(s2) 

+  a(l-a)5  (1-6)  4>e(ys1+8s2)  <J>e(s2) 

+  (1-a)  ( 1-ot)  5  ( 1-6 )  <J)E  (ys^)  q>E  (s2) 

+  aa(l-6)6c))E(s1+8s2)  ^(Bs^  4>e(YS2) 

+  (l-a)a(l-6)6<t>E(s1)  <D£(6s1)  4>e(ys2) 

+  a  ( 1-a)  (1-6)  6<}!e(s;l+3s2)  4>e(ys2) 

+  (1-a)  (1-a)  (1-6)  5(pE(s1)  1>e(ys2) 

+  aa  (1-6)  (1-5)  $e(ys1+3s2)  |>e(3s1)  ^£(ys2) 

+  (1-a)  a  (1-6)  (1-5)  $e(ys1)  ^(Bs-j.)  4>e(YS2) 

+  a  (1-a)  (1-6)  (1-6)  ^e(ysl  +  3s2)  1>e(ys2) 


+  (1-a)  (1-a)  (1-5)  (1-6)  ^{ys^  4>e(YS2) 


x  ( s ^ r s 2 )  =  f6$E(s2)  +  {l-6)  <J>e(ys2)  ] 

n+1 

x  [aa6$E(s1+6s2)  <?e(8s1)  +  (1-a)  a<5 <J>E (s1)  <j>E (Ss^) 

+  a  (l-a)  6  4>e(s1+Bs2)  +  (1-a)  (1-a)  S^s^ 

+  aa  (1-6  )  4>e(ys1+ds2)  <pe (Ss^) 

+  (1-a)  a  (1-6  )  $E  (YSj^)  $E  (Ss1)  +  a  (1-a)  (1-6)  (ys1+6s2) 

+  (1-a)  (1^)(1-6)4>e(ys1)  ] 


$ y  v  (s,,s9)  =  [64>r(s,)  +  (l-o)c{)Tr(YS0)  ]  (III. C. 7.1) 

Xn'  n+1  12  E  2  E  2 

x  [a$E (3s)  +  ( 1-a)  ]  x  ta6<j>E(s1+6s2)  +  (l-a)6(()E(s1) 

+  a  (1-6)  Jje(ys1+3s2)  +(l-a)  (1-6 )  4>E  (ys^)  ] 

For  the  special  cases  of  the  EMA(l)  process.  III. C. 7.1  reduces 
to  the  results  given  in  Lawrence  and  Lewis  [Ref.  5] . 

D.  THE  MOVING  MINIMUM  MODEL 
1 .  Introduction 

Another  possible  scheme  that  can  be  used  to  generate 
one-dependent  sequences  of  random  variables  with  marginal 
Exponential  distribution  is  the  so-called  minimum  model.  With 
this  model  the  {Xnl  sequence  is  generated  by  taking  the  moving 
minimum  value  of  two  Exponential  random  variables.  The  proposed 
generation  scheme  is 

=  MIN(En,bEn_1) ,  (III. D. 1.1) 
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where  {X^,  n  =  1,2,...}  is  a  sequence  of  random  variables 
with  marginal  Exponential  distribution,  (En,  n  =  0,1,...}  is 
an  iid  sequence  of  Exponential  random  variables  with  unit 
mean,  and  b  >_  0 .  This  will  produce  an  {X^}  sequence  with  a 
rate  of  and  an  expected  value  of  This  expected 

value  produces  one  difficulty  since  it  is  a  function  of  the 
parameter,  b.  This  complicates  comparisons  between  results 
with  different  parameter  values  and  decreases  the  value  of 
scatter  plots  and  sample  paths.  However,  this  difficulty 
can  be  easily  removed  by  multiplying  the  (xM  by  -g— .  The 
generation  scheme  then  becomes 

X  *  ^X'  =  MIN  (  [— ^— •]  E  ,  [b+1  ]  E  -.  )  ,  (III. D. 1.2) 

n  b  n  bn  n- 1 

with  (E  }  and  b  defined  as  before.  The  (X  }  has  a  rate  of 

n  n 

one  and,  hence,  an  expected  value  of  one.  This  facilitates 
comparisons  for  different  parameter  values  with  the  NEMA(l) 
discussed  in  III.C  which  produces  random  variables  with  unit 
means . 

The  investigation  of  the  moving  minimum  model  is  moti¬ 
vated  by  the  previous  result  in  III.C. 2  that  linear  additive 
models  have  a  constrained  range  of  serial  correlation.  The 
hope  is  that  the  non-linearity  of  the  moving  minimum  model 
will  obviate  this  constraint.  The  minimum  scheme  has  been 
used  by  Tavares  [Refs.  22  and  30]  to  generate  first-order 
autoregressive  exponential  processes  and  by  Marshall  and 
Olkin  [Ref.  31]  to  generate  correlated  bivariate  Exponential 


variables . 


1 


2 .  Correlation  Structure 

The  first-order  serial  correlation  can  be  computed 
by  the  following  approach 


E(Xn+lV 


E ( [MIN( [=g±]En+1, [b+l]En) ] 
[MIN (  [b+l]En_1)  ]  )  . 


The  terms  inside  the  expected  value  can  be  made  independent 
by  conditioning  on  the  value  of  En*  The  E(Xn+^Xn)  is  then 
found  by  multiplying  the  conditional  result  by  the  density 
of  En»  and  integrating.  Implementing  this  approach  we  have 


b+1 . 


E(Xn+lXn}  =  J  E ( [MIN ( I— p— 1 En+1 , [b+1] y ) ] 


b  n+1  ‘ 


[MIN  (  [— y  /  [b+HE^J  |En=  )e~ydy 


E<VlV  * 


(III. D. 2.1) 


rb  +  1 


x  ( E  [ MI N  (  [^]y,  [b+l]En_L)  ]  )e-ydy 


The  expected  value  of  the  minima  can  be  calculated  as  follows: 


h,,  (b+1)  y  .  -t ~ 

E(MIN[(^i)En+1,  (b+l)y])  =  /  x(bTT)e  dx 


00  bx 

+  /  (b+l)y  (mlyie  ^dx 

(b+1)  y  °  1 
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bx 

E(MIN[(^)En+1,(b+l)y])  =  -xe  ETI 


(b+1)  y 
0 


(b+1)  y  — 

+  /  e  dx  +  (b+1) ye  y 

0 


bx 


=  -  (b+1)  ye~ky  +  (^l)j'b+1)y(Ekl)e-b*l 

-by 


dx 


+  (b+1) ye 


E  (MIN  [  (~jj~)  En+1  /  (b+1)  y  3 )  =  ^(l-e“by) 


(III.D.2. 2) 


Similarly, 


E  (MIN [ (~£~) y , (b+1)  ER_ x ]  )  =  / 


(b+1)yyb  ,  l  -E7T. 


x(b+T)e  dx 


+  f  <b±i2X(  1  )e'b+1dx 

(b+1)  y/b  b  E*1 


-htt1 {b+1)y/b  <b+i)y/b  -E?r, 

=  -xe  !  +  J  e  dx 

|0  0 


(b+1) y  -y/b 
b 


,  a-y/b  +  (b+1)/<b+1)y/b(EiT)/B^dx 


.  (b^ily  -y/b 
b  e 
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(III. D. 2. 3) 


) 


E (MIN [  y  ,  (b+1)  ] ) 


=  (b+1)  (1-  e-y/b) 


Using  III. D. 2. 2  and  III. D. 2. 3  in  III. D. 2.1  produces 


E(Xn+lV 


=  /  (^A)  (l-e'by)  (b+1)  (l-e_y/b)e'ydy 


(b+1) 


b  +  1 
b 


(b+1>  +  — a±Hi 

b  ( 1+b+i) 


(b+1) 

b2+b+l 


Therefore,  since  E(Xn)  =  1 


COV(X  .  ,X  )  =  - 1 

n+1  n  b2+b+l 


b  +b+l 


and 


CORR(X„+1,Xn) 


b 

b2+b+l 


(III. D. 2. 4) 


Thus  the  model  allows  a  range  of  correlations  from 
[0,-j].  The  minimum  value  is  achieved  when  b  is  zero  or  in 
the  limit  as  b  tends  to  infinity.  The  maximum  value  is  achieved 
when  b  is  equal  to  one.  An  interesting  aspect  of  the  correla¬ 
tion  structure  of  the  moving  minimum  model  is  that  reciprocal 
values  of  b  produce  equal  correlations.  This  is  a  similar 
kind  of  "invertibility"  found  for  the  other  moving  average 
models  discussed  in  III.C.  The  range  of  b  could  be  restricted 
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to  the  interval  [0,1]  without  reducing  the  possible  range  of 
correlations.  However,  doing  so,  as  with  the  NEMA(l)  model, 
would  ignore  the  fact  that  in  the  non-normal  case  character¬ 
istics  other  than  correlation  may  be  different  in  the  case  of 
b  and  Also,  as  is  the  case  for  all  the  first-order  moving 
average  processes  addressed  in  this  paper,  the  correlation 
for  lags  greater  than  one  are  zero.  So  the  study  of  correla¬ 
tion  structure  is  limited  to  the  study  of  the  serial  correla¬ 
tion  with  lag  one. 

3 .  Negative  Correlation 

The  range  of  possible  correlations  can  be  extended  in 
a  fashion  similar  to  the  NEMA(l)  model  (see  III.C.2)  by  the 
use  of  correlated  or  antithetic  variables.  Using  this  approach 
the  generation  formula  becomes  for  antithetic  variables 


=  MIN(  [^]En,  [b+l]E^_i)  , 


(III. D. 3.1) 


where  all  variables  are  defined  as  in  III. D. 1.1  and  {E', 

n 

n  =  0,1,...}  is  generated  from  the  (E  }  sequence  using  the 

-E  n 

relationship  -  -  In  ( 1-e  n) .  Note  that  this  implies  that 
{ E^}  is  also  iid  Exponential  with  unit  mean.  Again 


E(VW 


E[  (MIN[  (££i)  En,  (b+l)E^_1]  ) 


x  (MIN ( (2£±) En+1 , (b+1) E^] ) ] 


and  cor.uitioning  on  the  value  of  E  ,  multiplying  by  the  density 
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) 


of  En,  and  integrating  produces 


E(XnXn+l)  =  /  E(  [MIN([^i]y,  [b+UE^)  1 


[min  {  En+1 ,  [b+1]  E^)  ]  |  En  =  y )  e  ydy 


-Y, 


E  {xnXn+l)  *  /  (E  [MIN  ([=£=•]  y,  [b+1]  E^)]) 


(III. D. 3. 2) 


(E [MIN ( En+1 , - [b+1] In [ l-e'y ] ) ] ) e'y  dy 


The  first  expected  value  is  identical  to  III. D. 2. 3.  Thus 


E(MIN[(~^)y,  (b+l)E^_1l) 


=  (b+1) (l-e~y/b) 


(III. D. 3. 3) 


The  second  can  be  calculated  as  before. 


E  (MIN  [  (~j~)  En+1,-(b+l)  In ( l-e~y)  ]  ) 


- (b+1) In ( l-e'y ) 


bx 

.  b  .  ~H+T. 

X^H+T  e  dx 


bx 

oo  -  !■ 

+  /  - (b+1) ln(l-e"y) e  D+idx 

-  (b+1) In ( l-e~y ) 


=  -xe 


bx 

'5+T 


-  (b+1)  ln( l-e'y)  - (b+1) ln(l-e'y) 


0  0 
-  (b+1) In ( l~e~y) ebln ( 1-6  Y) 


bx 

,  b  .  'b+1  . 
(b+T> e  dx 
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making  the  appropriate  changes  to  the  limits  of  integration 
produces 


E(VW 


Then 


'b+1)2-^-  (^i)2/  zbdz  + 


f°  u'<1+b)(l-u)b(-|i) 
b  u 


,b+l.  ,b+l*  b+1  ^  r1  1/b ,,  .b. 

("5_)  '  '  ~T  /  u  (1'u)  du 

D  0 

r (l+p) r (i+b) 

r  ( 2+b+g-) 

pr(^)br<b) 

(l+b+~)  (b+^)F{b+^) 

b2r(g) r (b) 

(b2+b+l)  (b2+I)  r  (b+^)  * 


C0VlXn'Vl>  " 


b2r  (g-)  r  (b) 
(b2+b+l) (b2+l) r (b+i) 

o 


-  1 


and 


C0RR'Xn'Vl>  " 


b2:(i)F(b) 
(b2+b+l)  (b2+l)  T  (b+J) 


-  1 


:iII.D.3.5) 


Like  the  expression  for  positive  correlation,  this 
expression  is  also  symmetric  with  respect  to  reciprocal  values 
of  the  parameter.  It  attains  a  minimum  value  of  minus  one- 
third  when  the  parameter  value  is  one.  Graphs  of  the 
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correlation  as  a  function  of  the  parameter  b  for  both  posi¬ 
tive  and  negative  correlations  are  provided  in  Figures  III. D. 3.1 
and  III. D. 3. 2,  respectively.  Unlike  the  NEMA(l)  model  which 
requires  an  additional  parameter  (see  equation  III.B.6)  to 
achieve  a  full  range  of  negative  correlations,  the  moving 
minimum  model  can  achieve  its  full  range  with  the  single 
parameter  b  and  an  antithetic  sequence  {E^}. 

4 .  Joint  Density  of  Xn  and  Xn | ^ 

Calculation  of  the  joint  density  of  Xfi  and  xn+^  is 
possible  using  a  conditioning  argument  to  determine 
P  ( X  <x|E  =  z)  and  P(X  .  ,<y|E  =  z)  .  These  values  along  with 
the  probability  that  En  takes  on  a  given  range  of  values  are 
sufficient  to  determine  the  joint  distribution  function  of  X 

n 

and  xn+2_*  The  form  of  the  distribution  will  vary  depending  on 
whether  one  is  above  or  below  the  line  X^+1  =  bXn.  The  joint 
density,  where  it  exists,  is  determined  by  differentiating 
the  distribution  function. 

From  III. D. 1.2  we  have 

Xn  =  MIN(  [^~]En,  [b  +  l]En_1)  . 


Then 
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P(z  £  REGION  2) 


P(z  £  REGION  3) 

P(z  £  REGION  3) 


Now  by  definition 


(III. D. 4. 4) 


(III. D. 4. 5) 


P(X  <X,X  <  V I z  £  REGION  i)  =  P(X  <x|z  e  REGION  i) 
n  —  n+i  —  1  n  — 

x  P  (X  i  <_  y  |  z  £  REGION  i)  (III.  C.  4. 6) 

because  when  conditioned  on  the  value  of  En,  these  probabili¬ 
ties  are  independent.  Using  the  above  equation,  III. D. 4.1, 

III. D. 4. 2,  and  the  definition  of  the  regions  in  Figure  III. D. 4.1, 


PlXnix-Xn+1^1 

Z  £  REGION 

1) 

=  1 

(III. D. 4 .7) 

X 

e«nix,Xn+liy: 

z  €  REGION 

2) 

.  b+1 

=  1  -  e 

(III. D. 4. 8) 

x  by 

P(Xnix'Xn+liy!z  4 

REGION  3) 

= 

, .  b+1 .  , .  b+1 . 
(1-e  ) (1-e  ) 

(III. D. 4 .9) 

Using  the  results  of  III. D. 4. 3  through  III. D. 4. 5  and  III. D. 4. 7 
through  III. D. 4. 9  we  can  compute  the  joint  distribution  of  Xn 
and  when  y  >  bx  by  using  the  relation 
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(III. D. 4 .10) 


P(Xn£x,Xn+1  <  y)  =  l  P(Xn  ix,Xn+1  <  y|  z  e  REGION  i) 


*  P (z  e  REGION  i) 


bx 


bx 


=  1  (l-e  b+l)+(l-e  b+1)  (e  b+1  -  e-^b+1 ) 


+  (l-e*^1)  e~y/ lb+^^ 


bx  bx 


_ _  _ (x+y)  y 

1-.  b+b+e_b+1-e~x-e  b+1+e  b+1  +e  ^ 


-y  -<ETT><x+tb+1>y) 
-  e  -e  y+e 


f«ni-.xn+1iy)  = 


,  -x  -v,  b+1 

l-e  -e  ■* +e 


(III .D. 


Similarly,  when  y  <  bx  (i.e.  when  you  are  below  the 
line  bx  =  y) ,  the  range  of  possible  z  values  can  be  separated 
into  three  regions.  See  Figure  III. D. 4. 2.  Then 


P(z  *£  REGION  1) 


P(z  £  REGION  2) 


P  (z 


(III.D. 


P(z  t  REGION  2) 


P(z  REGION  2) 


P(_Z_<  z  <  _£*_) 

^  b+1  2  -b+1'  ' 


(III.D. 


P(z  REGION  3; 


P  (Z  > 


bx  , 
b+1' 


4.11) 


4.12) 


4.13) 
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P(z  £  REGION  3) 


(III. D. 4 .14) 


Using  III. D. 4.1,  III. D. 4. 2,  III. D. 4. 6,  and  the  definitions  of 
the  regions  in  Figure  III. D. 4. 2,  the  following  results  hold 
for  the  given  region. 


P(X  <  xfX„.,  <y|z  REGION  1] 
n  —  n+1  —  -1  1 


(III. D. 4 .15) 


P  (X  <  X,  X  <  y  I  z  ?  REGION  2 )  = 

n  —  n+i  —  1 


=  1  -  e 


-bJL 

b+1 


(III. D. 4. 16) 


P(Xn < x,Xn+1  < y|z  £  REGION  3)  = 


(1-e  b+1) 


(III. D. 4 .17) 


x  (l-e'b+1) 


Combining  III. D. 4. 10  with  III. D. 4. 12  through  III. D. 4. 17  yields 
for  bx  <  y 


P(Xn<  x,Xn+1  <y)  =  l-e~y-e  x+e 


(III.D.4 .18) 


x 

.g-T-y 


1  -X  -V  b+T  1 

Let  F  (x, y)  =  1-e  -e  Y*e  ,  the  distribution 

function  of  X  and  X  ..  when  bx  <  v;  and  let  f^tx.y)  be  the 
n  n+i 

joint  density  of  Xn  and  XR+^  when  bx  <  y.  Then 


fl(x'Y)  =  hhFl{*,y) 


3  3  r,  -x  -y,  b+1  , 

=  —  ■*—  [1-e  -e  1  +e 

3x  oy 


3  ,  -y  S+T"y, 

-  ^[e  -e  1 
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f1(x,y) 


<bTTla 


X 

•s+ry 


bx  <  y?  x  >  0, 


(III.D.4 


2 

For  bx  ,<  y,  the  distribution  function,  F  '::,y),  is 
1-e  ^-e  x+e  b+1.  If  f2(x,y)  is  the  joint  density  of  Xn 
and  Xn+^  when  bx  >  y,  then 


f  (x,y)  = 


_ i (x  v)  *  - - r 

3x  3y  F  [X,y)  3x  3y  1  e  e  e 


-*-&r 


b+l. 


ax1 


b+l 


f  (x,y)  = 


-x--^L 

(s5r)e  b+1  y  <  bx:  sr  >  °. 


(III.D.4 


Note  that  there  is  a  positive  probability  that  the 
point  (xn'xn+]_)  lies  on  the  line  bx  =  y.  This  probability 
can  be  computed  as  follows.  We  have 


X 


rb+l 


n  '  MI“<(V1En'lb+11En-l> 


xn+l  “  MIN(  [^i]En+1,  tb+i:En) 


The 


b+l, 


point  (X  , X  lies  on  the  line  bx  =  y  when  X„  =  (-r— )E 

n  n+l  n  b  n 


and  X  ,.  =  (b+l)E  .  Now 
n+l  n 


p<xn-  'T'VVl-  Ib+1)En> 


'  p  <  ItT>  En  i  lb+1 1  En-1 !  lb+1 


.19) 


.20) 
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The  events  in  the  right  hand  side  can  be  made  independent  by 
conditioning  on  the  value  of  E  .  Then 


P(xn.  [TT]En;Xn+1=  [b*UEn> 


b2+b+l 


(III. D. 4 .21) 


Because  there  is  a  positive  probability  of  lying  on  the  line 
bx  =  y,  the  moving  minimum  model  can  be  said  to  have  a  line 
degeneracy.  An  important  implication  of  the  positive  proba¬ 
bility  of  (Xn,Xn+1)  lying  on  the  line  bx  =  y  is  that  the 
moving  minimum  model  will  produce  runs  of  values  of  constant 
ratio  b.  The  values  of  {Xn}  in  these  runs  will  be  decreasing, 
equal,  or  increasing  for  b  less  than,  equal  to,  or  greater 
than  one,  respectively.  The  length  of  the  runs  will  be  geo- 

metrically  distributed  with  parameter  — = -  for  the  positive 

b  +b+l 

correlation  case.  It  was  this  kind  of  degeneracy  in  the 
Exponential  autoregressive  model,  EAR ( 1 ) ,  that  proved  to  be 
one  of  the  model's  major  weaknesses.  The  degeneracy  also 
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occurs  in  the  Tavares  autoregressive  model  and  the  bivariate 
Exponential  pairs  derived  by  Marshall  and  Olkin. 


The  probability  of  lying  above  or  below  the  line  bx  =  y 
can  be  easily  found  by  integrating  the  appropriate  joint  den¬ 
sity  over  the  area  desired.  Thus,  for  bx  <  y. 
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5.  Conditional  Expectation  and  P (X  >  X  ) 

- E -  n  +  1  n 

Besides  the  correlation  coefficient,  there  are  two 

other  chaiacterizations  of  the  joint  distribution  of  Xn  and 

Xn+1  which  we  have  considered.  They  are  the  conditional 

expectations  and  the  P(Xn+^  >  xn)  •  Both  of  these  quantities 

can  be  derived  by  considering  the  four  possible  sets  of  values 

for  Xn  and  X^+^ ,  computing  the  probability  of  each  set  occurring, 

and  weighting  the  conditional  expectation  or  probability  by 

its  probability  of  occurrence. 

First,  the  probability  of  occurrence  for  each  set  of 

values  must  be  calculated.  Consider  the  case  where  X  =  (^rr^)E 

n  on 

and  Xn+1  “  (~b~i  Ea+i ‘ 


p<Xn*  ^".'VrfT'vi1 


*  En-1  •  'TT1  En+1  i'W]V 


By  conditioning  on  the  value  of  En,  the  events  on  the  RHS 
become  independent.  The  calculation  then  proceeds  in  a 
straightforward  way . 


P,Xn=  <TT1En'xn+l-  >T'Vi: 


/oP([^]yi  [b+l]En_x,  [2±i]En+1  <  [b+l]y!En  =  y)e'ydy 


^P(En-l>b)P(En+l-by)e'ydy 
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inspection. 
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Case 


X 


X 


X 


X 


n 


n 


n 


n 


iT'V^r  't'Vi 


Ib+1'En 


b+l. 


lb+llEn-l'xn+l*  ^Vl 


[b+1]En-l'Xn+l  =  fb+l]En-l 


b+l 

b 

by 

b+l 

b 

b+l 


Weighting  these  conditional  expectations  by  the  probabilities 
in  III. D. 5.1  through  III. D. 5.4  yields  the  final  result. 


E(Vilxn-r»  * 


i  E(X  .  | X  = y ;case  i)P(case  i) 
i=l  n  x  n 


-  [• 


(b+l) (b  +b+l) 


? - ]  +  (by)  (-^ 

‘  it  il  \  1 _  ^ 


b  +b+l 


+  (^~)  (-5^ - )  +  (b+l)  (- 

"  b  +b+l 


[b+l] [b  +b+l ] 


-) 


E(Xn+liXn  =  ^  = 


i  + 

b  +b+l 


(III. D. 5. 5) 


It  is  quite  surprising  that  the  regression  of  Xn+^  on  Xn  is 
linear  in  y,  considering  the  non-linearity  of  the  process  which 
generates  the  {Xn) . 

The  conditional  expectation  of  X^  given  Xn+1  can  be 
derived  with  equal  dispatch. 
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Case 


£<xnlxn+1ii± 


x„  =  lT'8n'Vl=  'T'Vl 


Xn*  <— 1En'Vl=  lb+11En 


Xn=  tb+1 1  En_i  >  xn+i  = 


Xn  =  <b+1>En-l'X„+i  =  Ib+1JEn 


Using  III. D. 5.1  through  III. D. 5. 4  as  before 


E(XniXn+l  =  ^  = 


1  E(X  |Xn+i  = y,case  i)P(case  i) 
i=l 


(b±l)  ( - ^ - )  +  (Z)  (—^ - ) 

D  [b+1] [bz+b+l]  °  o  +b+l 


+  (b+1)  (-=^ - )  +  (b+1)  ( - — - 

b^+b+l  [b+1] [d  +b+l ] 


E(Xn!Xn+l  =  ^  " 


(III. D. 5. 6) 


b  +b+l 


The  probability  that  Xn+^  is  greater  than  Xn  can  also 
be  easily  computed  if  one  is  careful  to  differentiate  between 
the  case  where  b  <  1  and  b  >  1 . 


Case 


^n+l^V 


Xn=l^VXn+l=[^En+l 
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1 


Xn"lVIEn'Vl*[b+11En 

Xn  =  lb“'Vl'Xn*l-|bT|E,H 

x„  =  <b+11En-l'Xn+l=  Ib+11En 
Thus  when  b  <  1  we  have,  using  III. D. 5.1  through  III. D. 5. 4 


0  if  b  <  1, 
1  if  b  >  1 


b+1 


1 

2 


P(X„  >  X) 
n+l  n 


i=l 


P  (X  , ,  >  X  lease  i)P(case  i) 
n+i  n 


[b+1] [bz+b+l] 


-)  +  (err) 


b+1'  V+b+l 


+  (|)  ( 


[b+1] [b  +b+l] 


-) 


P(Xn+l>Xn) 


1 

2 


(b+1) (b^+b+l) 


•,  b  <  1, 


(III. D. 5. 7) 


A  similar  computation  with  b  >  1  again  using  III. D. 5.1 
through  III. D. 5. 4  yields 


f'Vi’V 


(b+1) (b  +b+l) 


b  >  1 


(III. D. 5. 3) 


Thus  a  graph  of  P(Xn+1  >  X  )  will  have  a  discontinuity  at  b  =  1 
when  case  2  switches  from  a  probability  of  zero  to  a  probability 
of  one.  This  graph  is  presented  as  Figure  III. D. 5.1.  The 
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minimum  value  of  one-third  occurs  at  b  =  1.  The  maximum  value 
of  two-thirds  occurs  at  b  =  1+.  The  moving  minimum  model, 
therefore,  has  a  greater  range  of  values  for  the  P(Xn+^>Xn) 
than  does  the  NEMA(l)  model.  However,  the  greater  range  for 
the  p(xn+i>xn)  and  greater  range  of  correlations  must  be 
balanced  against  the  degeneracy  of  the  model. 

As  was  noted  with  the  NEMA(l)  model,  the  correlation 
in  non-normal  models  does  not  define  the  joint  properties  of 
Xn  and  xn+i*  Although  the  cases  of  b  and  g  are  indistinguish¬ 
able  from  the  viewpoint  of  correlation  (see  III. D. 2. 4  and 
III. D. 3. 5),  these  cases  will  have  significantly  different 
sample  paths  as  indicated  by  III. D. 5. 7,  III. D. 5. 8,  and  the 
discussion  of  runs  up  and  down  in  III.D.4. 

Three  examples  of  sample  paths  for  different  b  values 
are  given  in  Figures  III. D. 5.2  through  III. D. 5. 4.  The  degen¬ 
eracy  of  the  model  is  clearly  present  in  the  sample  paths  as 
a  tendency  to  produce  runs  of  equal,  increasing  or  decreasing 
values,  respectively.  A  comparison  of  Figures  III. D. 5. 3  and 
III. D. 5.4  quickly  demonstrates  that  while  these  two  sample 
paths  have  the  same  correlation,  they  produce  significantly 
different  {Xn}  sequences.  This  is  a  graphic  indication  that 
non-normal  processes  are  not  determined  solely  by  their 
correlation  structure. 

Figures  III. D. 5. 5  through  III. D. 5. 7  are  the  scatter 
plots  associated  with  the  sample  paths  in  Figures  III. D. 5. 2 
through  III. D. 5. 4,  respectively.  Here,  too,  the  degeneracy 
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of  the  model  is  clearly  present  in  the  tendency  of  the 
(Xn,Xn+i>  plots  to  lie  on  the  line  Xn+1  =  bXn>  The  slope 
of  this  line  determines  whether  the  runs  are  of  equal,  in¬ 
creasing  or  decreasing  value. 

6.  Conditional  Expectation  and  P (X  .  >  X  )  for 
Antithetic  Variables'  n  n 

Results  similar  to  those  obtained  in  III.D.5  can  be 
obtained  for  the  moving  minimum  model  with  negative  correlation. 
The  procedure  for  determining  the  conditional  expectations  and 
the  probability  that  Xn+^  is  greater  than  Xn  using  antithetic 
variables  is  exactly  the  same  as  that  in  the  previous  section. 
First,  the  probability  of  each  of  the  four  possible  combina¬ 
tions  of  Xn  and  values  is  computed,  the  conditional 

expectation  or  probability  is  computed  for  each  case,  and 
the  final  weighted  sum  of  conditional  expectations  or  proba¬ 
bilities  is  finally  computed.  In  one  instance  no  closed  form 
answer  is  available  and  numerical  procedures  are  used. 

Recall  that  the  generation  scheme  when  using  antithetic 
variables  is 

X„  =  MIN(  [^]E.  [b+l]E’  ]  ,  (III.  D.  6.1) 

n  on  n— l 

where  (En,  n  =  0,1,...}  is  an  iid  sequence  of  Exponentially 

distributed  random  variables  with  unit  mean,  {E',  n  =  0,1,...} 

n 

is  generated  from  the  {ER}  sequence  by  the  relationship 
“  E 

E^  =  -  In  (1-e  n)  which  implies  that  {En}  is  also  iid  Exponen¬ 
tial  with  unit  mean,  b  >  0. 
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First,  consider  the  case  where  X  =  and 

n  on 

Vi  -  'T'Vr  Then 


P(x„-  (^IEn,xn+1=  [Sii]En+1i 


-  p'‘6Ei)Eni  'b+1)EA-i'  'x’Vii  (b+1)EA> 


Using  the  standard  conditioning  argument  produces 


P(Xn*  t65iIEn.Xn+1=  l^i]En+1) 


=  J  p(t^]yi  [b+l ] E^,  [^i)En+1  <  -[b+l]ln(l-e"y)  |En  =  y)e‘ydy 


•  «  rb-Hi 


-  /  P  (E^  >  J)  P  (En+1  <  -bln  [l-e*y] )  e“ydy 


/  e-y/b(l-[l-e-y]b)e-ydy 
0 


-  -(i+§) y 


/  e 
0 


oo  (1+—) 

dy  -  /  (1-e  y)b(e  y)  b  dy 
0 


The  first  integral  is  straightforward.  In  the  second  integral, 
the  change  of  variable  u  -  e  y,  =  dy  makes  this  integral 

recognizable  as  the  integral  of  a  Beta  random  variable.  Making 
this  change  of  variable  and  making  the  appropriate  changes  in 
the  limits  of  integration  produces 
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P(Xn=  [TT1En-Xn+1=  [b+liVl> 


=  bTT  "  /  d-u)b(u)bdu 


b 

b+1 


r (i+b)r (i+£) 
r  (2+b+^) 


Plx„=  iriVVi'^'W 


b 

b+1 


b2r (b) r  (|) 


(b2+b+l) (b2+l)r (b+i) 


(III. D. 6. 2) 


In  the  second  case,  X  =  [^r^]E  and  X 

n  b  n  n+1 

ceeding  as  before 


[b+1 ] . 


Pro- 


P(xn-  1*^1  VVi*  tbH-1  J  E^> 


-  P“TT-lE„i  [b*11EA-l'  Ib+1lEAi  1TTlEn+l) 


/oP([^]y  <  [b+l]Ej^_1,-[b+l]ln[l-e~y]  <  En+]_  | Efl  =  y)  e_ydy 


/  P (EA-1  >  b} P (En+l  >  ~bln [1~e~Y] > e_ydy 


/  e-y/b(l-e-y)be-ydy 
0 


00  h  -<1+h)y 

/  (1-e  y)be  dy 

0 
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This  is  the  same  integral  as  the  second  integral  in  the  first 
case.  Thus 


,b+l 


P(Xn=  (V^n-Vr 


b2r(b)  r(|) 


(b2+b+l)  (b2+l)  r  (b+i) 


(III.  D.  6. 3) 


Next  consider  the  case  where  X  =  (b+l)E‘  .  and 

n  n-1 

Vi  -  'x'Vr  Then 


P(xn-  [Miis;.rxntl-i^]Entl) 


•  E<  »+1)EAi  V  i  lb+11EA 


/  P([b+1]E^<  t^]y,  C^lEn+ii  fb+l]ln[l-e“y]  |En  =  y)e~ydy 


/  p(En  Ib}  P(En+l  -  -bln[l-e"y] )  e-ydy 


/  e  ydy  -  /  e  ^  ey-/ey  (-e  y)  ^dy  +  /  (1-e  y)  ^  (e  y)  ^ 
0  0 


-d+R)y 

« 

0 


dy 


The  first  two  integrals  do  not  present  a  problem.  Making  the 
change  of  variable  z  =  l-e~y,  dz  =  e  ydy  makes  the  third  inte¬ 
gral  easy.  The  last  integral  is  the  same  as  the  second  integral 
in  the  first  case.  So 
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p<xn=  [b+l]E’_i,xn+iMb+l]E-) 


b2r(b)r(|) 

(b2+b+l)  (b2+l)  r  (b+jj) 


(III. D. 6. 5) 


The  conditional  expectations  given  a  specific  case 
of  the  values  of  Xn  and  Xn+1  can  be  written  by  inspection. 
Hence, 


Case 


S«n+1lXnJLii 


xn’  'TIE.'Vr  'TT!W 


b+1 
b  * 


xn-  irIE.'Vr  tb+llEA! 


-  (b+1) In (1-e  b+1) 


Xn=  tb+1]EA-l'Xn+l  =  [VlEn+r 


Xn=  [b+UBA-rVi-  Ib+1]EA; 


b+1. 


Combining  these  results  with  equations  III. D. 6. 2  through 
III.D.6.5  and  letting 


G  = 


b2r(b)  r(jj) 

(b2+b+l) (b2+l) r (b+jj) 


we  have 
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E(xn+1lxn  =  y)  =  i  E(x„+1  Ix^-y.case  i)P(case  i)  (III. D. 6. 6) 


i=l 


n+1 1  n 


->Y 

(glj-  G)  -  G(b+l)ln(l-e  b+1) 


E(Xn+llXn=sy)  - 


+  G(^)  +  (bTT“  G)  (b+1) 


Jbsl 

2  -  G (b+1) (1+ln [1-e  b+1] ) 


(III. D. 6. 7) 


Similarly,  we  can  derive  the  expression  for  E(xn|Xn+1  =  y) 


Case 


Xn=[6Bi)En-X„+l=  'T'V!' 


Ia=[T1In'VrME;: 


b+1 
b  ’ 


I^UnU-e"5*1) 


xn*  lwllE;-i'Vr  if'Vr 


(b+1) 


Xn=  fb+1]EA-l'Xn+l=  [b+ljEA; 


(b+1) 


Then  using  III. D. 6.2  through  III. D. 6.5  and  again  letting 
G  = 


b2r (b) r (i) 


(b2+b+l) (b2+l)r(b+g) 


E  (xn  I  xn+i  =  y)  =  1  E(Xn|Xn4-1,case  i)  P  (case  i) 


i=l 


n 1 ~n+l 1 


__Z_ 

=  (^)  (g|r-G)-G(^ln(l-e  b+1)  +  (b+1)  G+  (b+1)  (g^-G) 
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-y0  b  -yo 

Hence,  if  we  can  find  yQ  such  that  (1-e  )  =  e  ,  then 

-y0 

P (Xn+1 >  xn)  =  1-e  .  The  required  solution  can  be  found  by 

numerical  means  for  any  given  value  of  b.  A  computer  program 
to  determine  y^  to  an  accuracy  of  10  ®  for  a  given  value  of 
b  and  to  compute  P(Xn+^  >  xn)  was  prepared.  A  graph  of  the 
results  is  presented  as  Figure  III. D. 6.1.  When  using  anti¬ 
thetic  variables,  the  moving  minimum  model  has  a  restricted 
range  of  possible  values  for  P(Xn+^>Xn).  The  maximum  value 
of  approximately  0.509  occurs  at  about  0.30.  The  minimum 
value  of  approximately  0.491  occurs  at  about  3.33. 

This  small  range  of  values  for  the  P(Xn+^>Xn)  Is 
shown  in  the  relative  indistinguishability  among  the  sample 
paths  displayed  in  Figures  III.D.6.2  through  III. D. 6. 4.  Of 
more  interest  are  the  scatter  plots  presented  in  Figures  III. D. 6. 5 


through  III. D. 6. 7.  In  these  plots  the  degeneracy  of  the  moving 
minimum  model  is  clearly  displayed.  When  Xn  achieves  a  value 
of  y  based  on  E  ,  then  X  . ,  is  constrained  to  have  a  value  less 
than  -ln(l-e  n) .  In  the  case  where  equality  is  achieved,  the 


second  case  in  the  discussion  in  this  section,  the  point  (X  , 

-X  ~x  i  n 

Xn+^)  lies  on  the  curve  e  n  +  e  =1.  This  constraint 

is  clearly  evident  in  the  scatter  plots.  Thus  the  moving 


minimum  model  displays  a  degenerate  behavior  for  negative 
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E.  THE  BETA- EXPONENTIAL  MODEL 
1 .  Introduction 

A  third  method  that  can  be  used  to  generate  correlated, 
marginally  Exponentially  distributed  random  variables  is  a 
special  case  of  the  Beta-Gamma  model  given  in  Lawrance  and 
Lewis  [Ref.  18].  This  model  generates  an  {Xn}  sequence  using 
the  relation 


xn  ■  Bn(q,1_q)En  +  Bn (1“q,q) En-1  n  =  1'2'***' 


(III. E. 1.1) 


where  (Bn(q,l-q),  n  =  1,2,...}  is  an  iid  sequence  of  Beta 
random  variables,  (Bn(l-q,q),  n  =  1,2,...}  is  an  iid  sequence 
of  Beta  random  variables,  (En,  n  =  0,1,...}  is  an  iid  sequence 
of  Exponential  random  variables  with  unit  mean,  {Bn(q,l-q)}, 
(Bn(l-q,q)},  and  {En}  are  mutually  independent,  and  0  <  q  <  1. 
The  density  for  a  Beta(m,n)  variable  is 


T (m+n) 
r  (m)  f (n) 


m-1.,  .  n-1 

x  (1-x) 


0<x<l;m>0;n>0. 


(III. E. 1.2) 


In  practice  the  Beta  random  variables  can  be  generated  from 
two  Gamma  distributed  random  variables  using  the  relationship 
B(m,n)  =  x r- rM- t  from  Xotz  and  Johnson  [Ref.  19],  where  G(K) 

Vj  ial )  \  n ) 

is  a  Gamma  random  variable  with  a  shape  parameter  of  K  and 
a  scale  parameter  of  one. 

This  is  a  special  case  of  the  Gamma  model  considered 
in  Chapter  II  of  this  thesis.  It  works  because, as  described 
by  Lewis  [Ref.  10],  in  III. E. 1.1  the  product  of  the  Bn(q,l-q) 
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variable  and  the  En  variable  is  a  Gamma  (q)  variable.  Simi¬ 
larly,  the  product  of  the  Bn(l-q,q)  and  the  En_^  variables  is 
a  Gamma  (1-q) ,  independent  of  the  Gamma  (q)  variable.  Conse¬ 
quently,  their  sum  is  an  Exponential  variable,  Xft.  The  (Xn) 
process  is  clearly  one-dependent,  as  for  the  NEMA(l)  process. 

Because  of  a  lack  of  a  closed  form  solution  for  the 
integral  of  the  Beta  density  when  the  limits  of  integration 
are  not  from  zero  to  one,  this  model  is  the  least  tractable 
of  those  considered  in  Chapter  III  of  this  thesis.  However, 
its  correlation  structure  can  be  determined,  an  expression 
for  the  Laplace-Stielt jes  transform  of  a  sum  of  r  random  varia¬ 
bles  can  be  derived,  and  the  probability  of  Xn+^  being  greater 
than  Xn  can  be  simulated.  An  advantage  of  this  model  is  that 
it  extends  directly  to  moving  average  Gamma  processes  (see 
Chapter  II) .  This  extension  is  not  possible  with  the  NEMA(l) 
or  the  moving  minimum  model . 

2 .  Correlation  Structure,  Positive  and  Negative 

The  correlation  structure  can  be  determined  using  a 
standard  approach.  We  have  using  III. E. 1.1 

XnXn+l  =  CBn(q,1_q)En+Bn(1_q'q)En-l) (Bn+l(q,1_q)En+l 
+  Bn+l(1“q'q)En) 

=  Vl(q'1-q)Bn(q'1-q)En+lEn+Bn+l(q'1'q)Bn(1-q'q)En+lEn-l 
+  Bn+]_  ( 1-q ,  q)  Bn  ( q  ,  1-q )  En+Bn+1  ( 1-q , q) Bn ( 1-q , q) 
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Taking  expectations  and  using  the  iid  nature  and  independence 
of  the  {Bn(q,l-q)},  (Bn(l-q,q)},  and  { En }  sequences  yields 

E(XnXn+l)  =  +  qd-q)  +  2q(l-q)  +  (1-q)2 

Hence, 

C0v(xn'xn+i)  = 


and 


CORR(Xn,Xn+1)  =  q (1-q) ,  0  <  q  <  1.  (III. E. 2.1) 

As  with  the  other  linear  additive  models,  this  correlation  is 
double  valued  and  lies  in  the  range  (0,^) . 

The  range  of  possible  correlations  can  be  extended 
to  negative  values  by  modifying  the  generation  formula  by 
including  an  {E^}  sequence.  Thus 

Xn  =  Bn(q'1_q,En  +  Bn(1"q,q) En-1'  (III. E. 2. 2) 

where  all  random  variables  and  constants  are  as  defined  for 
III. E. 1.1  and  { E n  =  0,1,...}  is  an  iid  sequence  having  a 
specified  correlation  with  the  {En}  sequence.  In  particular, 

E  and  E'  may  be  an  antithetic  pair.  The  correlation  of  the 
(Xn)  using  II.E.2.2  can  be  determined  in  the  same  way  as  before. 
Consequently, 
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(Bn(q,l-q)En  +  Bn(l-q,q)E^_1) 


XX  , , 
n  n+1 


(Bn+i(q'1_c3)En-n  +  B^i  (i-q/q)E;> 


n+l  "  n+1 


=  Sn+1  (q,  1-q)  Bn  (q,  1-q)  En+1En  +  Bn+1  <q,  1-q)  Bn (l-q,q)  E^E^ 


+  Bn+l(1-q'q)Bn(q'1-q)EnEA  +  Bn+l{1-q'q)Bn(1-q'q)EAEA-l 


Taking  expectations  as  before 


E(XnXn+l> 


q2+q<l-q)+q(l-q) [COV(En,E^) +l]+(l-q) 2 


1  +  q(l-q)COV(En,E^) 


Therefore, 


COV{Xn,Xn+1)  =  q(l-q)COV(En,Ei;) 


and 

CORR(Xn,Xn+1)  =  q(l-q)CORR(En,E^) ,  0  <  q  <  1.  (III. E. 2. 3) 

When  E'  =  E  ,  the  correlation  is  one  and  III. E. 2. 3 
n  n 

reduces  to  III. E. 2. 2.  When  E„  and  E'  are  antithetic  the 

n  n 

correlation  is  -0.6449  and  negative  correlations  result.  When 

q  is  0.50,  the  correlation  for  the  {Xn>  sequence  falls  in  the 

(-0.16,0.25)  range  depending  on  the  correlation  between  En 

and  E'. 
n 
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3 .  Laplace-Stielt j es  Transform  of  a  Sum 


t  =  y  x. , 

r  i=i  1 


(III. 3.1) 


where  {X^}  are  defined  by  III. E. 1.1.  Then 


Tr  =  X1  +  X2  +  *  *  *  +  Xr 


=  Bx (q,l-q)E1+B1(l-q,q)E0+B2 (q,l-q)E2+B2 (l-q,q)E1 


..  +  Br(q,l-q)Er+Br(l-q,q)Er_1 


=  Br(q, 1-q) Er+B1(l-q,q)E0+[B1(q,l-q)+B2 (l-q,q) ]E1 


t  t  ®  ^  ( g[  /  i—<3 )  j-  ^  ^  — '  *5 )  J  E^_  ^ 


-sT 


Let  $T  =  E(e  r) .  Then  using  the  iid  nature  and  independence 
r 

of  {Bn (q, 1-q)  } ,  {Br(l-q,q},  and  { En > 


-s  [Br  (q,l-q)  Er+B1  (l-q,q)  (q,l-q)  +B2(l-q,q)  }E1 

+...+{B  . (q, 1-q) +B  (l-q,q) }E  .] 

*T  (s)  =  E(e  r_i  r  r  1 

~sB  (q, l-c) E  -sB..  ( l-q,q)  E 

-  E(e  r  n)E(e  1  °) 


■s  t (B,  (q, 1-q) +B_ (l-q,q]E. 
x  E(e  1  *  1) 


-s[B  , (q,l-q)+B  (l-q,q) ]E  , 
x  . . .E(e  r  1  r  r  X) 
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4>m  (S) 
r 


(ir?iq[r7i)1'q<’,,BE(s>)r'1 


*  ‘IT?1  I'fBE(,,ir'1' 


where 

-s [B. (q,l-q) +B .  ,(l-q,q)]E, 

^(s)  =  E(e  1  11  1-1  )  . 

The  Laplace  transform  of  the  sum  of  two  Beta  random  variables 
is  a  confluent  hypergeometric  function.  Its  form  is  too  compli 
cated  to  be  of  significant  value  in  deriving  the  analytic 
behavior  of  the  Beta-Exponential  model . 

4.  Empirical  P(Xn|1  > Xn) 

Because  of  the  presence  of  the  Beta  random  variables, 
the  probability  of  Xn+^  being  greater  can  not  be  analytically 
determined  with  a  reasonable  amount  of  effort.  In  an  attempt 
to  establish  a  range  for  this  probability,  a  simulation  was 
used.  In  order  to  achieve  a  precision  of  at  least  0.001, 
sixty-eight  thousand  comparisons  were  generated  for  each  of 
ten  random  number  seeds.  The  Beta  random  variables  were 
generated  using  the  Kotz  and  Johnson  [Ref.  32]  relation 

f  —  \ 

B(m,n)  =  (n)"  explained  in  III.E.l.  The  Exponential 

sequences  were  generated  by  a  call  to  a  standard  generator  of 
Exponentials.  When  the  simulation  was  run  for  nineteen  values 
of  q  from  0.05  to  0.95  in  steps  of  0.05,  the  P(Xn+^>Xn)  was 
0.500  for  all  values  of  q. 
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Although  the  empirical  probability  that  Xn+^  is 

greater  than  X  is  constant  at  a  value  of  0.500,  reminiscent 
n 

of  the  autoregressive  model  of  Chapter  II. B. 6,  the  distribu¬ 
tion  of  Xn+^  -X  is  not  symmetric  and  no  simple  proof  for 
this  situation  has  been  found. 

The  low  serial  correlation  and  the  apparent  invaria¬ 
bility  of  the  p(xn+]_>xn)  stakes  the  use  of  sample  paths  and 
scatter  plots  of  little  value.  Samples  are  provided  in  Figures 
III. E. 4,1  through  III. E. 4. 12. 
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FIGURE  III 
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IV.  PRELIMINARY  DATA  ANALYSIS 

A.  INTRODUCTION 

During  the  period  1955  through  1969  a  weather  ship  sta¬ 
tioned  in  the  Gulf  of  Alaska  (50°N,145°W)  collected,  among 
other  data,  wind  speed  data  at  three  hour  intervals  [Ref.  33]  . 
The  existence  of  the  wind  speed  data  was  brought  to  Professor 
Lewis '  attention  when  a  student  in  Oceanography  asked  him  to 
provide  a  model  suitable  for  simulating  wind  velocity  data. 

The  simulated  data  was  required  as  input  to  models  of  ocean 
temperature  mixing.  A  copy  of  fifteen  years  of  wind  speed 
data  was  obtained  for  this  thesis.  The  intent  was  to  do  a 
preliminary  data  analysis  and  then  determine  whether  any  of 
the  models  discussed  in  this  thesis  could  provide  an  adequate 
representation  of  this  data  and,  hence,  a  method  for  gener¬ 
ating  wind  velocity  sample  paths  for  oceanography  simulation. 

Initially,  the  models  discussed  here  are  strong  a  priori 
candidates  for  data  of  this  nature.  Intuitively,  there  is  a 
strong  feeling  that  an  assumption  of  independence  among  the 
data  is  unwarranted.  Hence,  autoregressive  and  moving  aver¬ 
age  models  which  are  designed  to  reflect  the  behavior  of 
correlated  data  should  be  considered  likely  candidates.  The 
non-negative  nature  of  the  data  mitigates  against  the  use  of 
classical  time  series  anlaysis  which  is  based  on  the  assump¬ 
tion  of  a  normal  distribution,  and  hence  negative  values,  for 
the  series.  The  existence  of  zeros  in  the  data  tends  to  make 
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the  use  of  transformations,  like  taking  the  log,  somewhat 
less  appealing  than  otherwise.  These  considerations  indi¬ 
cate  that  the  models  discussed  in  this  thesis  should  be 
considered  as  likely  candidates  for  modeling  the  wind  speed 
data. 

B.  ANALYSIS  OF  THE  RAW  DATA 

There  were  43,800  data  points  available  for  analysis, 

2920  for  each  of  the  fifteen  years  between  1955  and  1969 
inclusive  (the  extra  data  for  leap  years  was  discarded) . 

Since  this  size  data  base  made  it  inconvenient,  if  not 
impossible,  to  manipulate  by  hand,  each  year's  data  was 
plotted  as  a  means  to  promote  familiarity  with  the  data. 

The  plot  of  each  year's  data  and  the  plot  of  the  data  averaged 
over  all  fifteen  years  (e.g.,  all  data  taken  at  0300  on  1 
January  of  each  year  were  averaged)  are  presented  in  Figures 
IV.B.la  through  IV.B.lp.  Several  characteristics  can  be 
observed  from  these  figures.  First,  and  perhaps  most  obvious, 
is  the  expected  yearly  cycle  of  the  data.  Values  at  the 
beginning  and  end  of  the  year  tend  to  be  higher  than  those 
in  the  middle.  Second,  the  data  is  discretized  to  a  large 
extent.  There  exist  obvious  horizontal  lines  of  equal  valued 
data.  A  sort  and  plot  of  the  entire  data  set  reveals  that 
the  data  consists  of  values  that  are  integral  multiples  of 
1.03  with  a  few  values  between  these  multiples.  Next,  on 
some  occasions  a  series  of  high  values  will  all  be  equal, 
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indicating  that  some  clipping  may  have  occurred  (see  Figure 
IV.B.la  vicinity  of  2575  and  2625,  Figure  IV.B.lb  vicinity  of 
400  and  525,  etc.).  These  last  two  characteristics  indicate 
that  statistical  properties  which  are  sensitive  to  the  be¬ 
havior  of  the  "tail"  of  a  distribution  may  be  affected.  The 
final  observation  about  the  data  is  that  there  are  apparently 
intervals  when  the  data  was  not  actually  collected.  These 
instances  appear  as  reasonably  long  strings  of  values  which 
have  a  strong  linear  appearance  (as  though  the  values  were 
produced  by  linearly  interpolating  between  two  boundary 
values) .  See  Figures  IV.B.lh  (vicinity  2400) ,  IV.B.lj 
(vicinity  2250),  IV.B.lk  (vicinity  50  and  1750),  and  IV.B.lm 
(vicinity  150  and  2725) . 

The  cyclical  nature  of  the  data  is  somewhat  more  apparent 
in  the  plot  of  the  data  averaged  over  the  fifteen  years  (see 
Figure  IV.B.lp) .  Additional  evidence  of  this  yearly  cycle 
is  presented  in  Figure  IV. B. 2.  This  figure  presents  twelve 
box  plots,  one  for  each  month.  The  data  values  plotted  are 
the  monthly  average  wind  speed  for  each  of  the  fifteen  years. 
The  interquartile  range  and  extreme  values  are  shown  in  a 
standard  fashion.  As  an  adjunct  to  this  analysis  of  the 
year  cycle,  the  coefficient  of  variation  for  the  monthly 
averages  was  computed.  The  coefficient  of  variation  was 
essentially  constant.  See  Table  IV.B.l.  This  will  have 
an  impact  on  the  choice  of  the  type  of  model  used  to  model 
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Box  plot  for  average  wind  speech  for  each  month  for  15  years 
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This  yearly  cycle  is  also  shown  in  the  periodogram  and  the 

log  of  the  periodogram  of  the  data  (averaged  over  15  years)  as 

presented  in  Figures  IV. B. 3  and  IV. B. 4,  respectively.  The 

periodogram  is  computed  from  the  data  in  the  following  way. 

N 

Let  { X  ,  n  =  1, 2 , . . . ,N}  be  the  raw  data  and  let  X  =  J  X.  be 

2  i=l  1 

the  mean  of  the  {X  }  sequence  and  at  be  the  variance.  Let  the 

n  x 

{ Y^ ,  n  =  1,2,...,N}  be  formed  from  the  {Xn)  sequence  using  the 
relation 


Y  =  X  -  X,  (IV.B.l) 

n  n 

where  N  =  2920  is  an  even  number.  The  Fourier  transform  of  the 

{Yn)  sequence  will  have  both  a  real  and  complex  component  and 

will  have  j  elements.  Let  {Zn,  n  =  l,2,...,j}  be  the  Fourier 

transform  of  the  {Y„}  sequence  and  let  Z._  and  Z.T  be  the  real 

n  3  ■K  j x- 

and  imaginary  components  of  the  jth  element  of  (Znl,  respectively. 

h 

Let  Pj  be  the  j  element  of  the  periodogram.  Then 

Pj  =  (Z^R  +  ZjI)/27TNo^,  j  =  1,2, ...,|  (IV. B. 2) 

defines  the  periodogram  of  the  {Xn}  sequence. 

The  periodogram  dramatically  presents  the  yearly  cycle 
( j  =  1)  as  the  dominant  effect  (P^  >  150) ,  although  there  is 
some  indication  of  a  six  month  cycle  (j  =  2,  P^  ~  9.0)  .  Some¬ 
what  surprising  is  the  apparent  lack  of  any  strong  time  of  day 
effect.  The  log  periodogram  reinforces  the  dominant  role  of 
the  yearly  cycle  and  indicates  that  six  month  and  six  and 
twelve  hour  cycles  (j  =  2,  j  =  1460,  j  =  2920  respectively) 
may  be  important. 
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second  harmonic  at  w.  (six  months). 


The  correlation  structure  of  the  data  is  presented  in  Table 
IV. B. 2.  The  column  indicated  as  "15  Yr  Avg"  is  the  average  of 
the  values  of  the  fifteen  years.  The  column  indicated  as  "15 
Yr  SD"  provides  the  standard  deviation  of  the  values  about 
their  mean.  It  is  not  the  standard  deviation  of  the  average. 
This  latter  quantity  can  be  obtained  by  dividing  by  the  square 
root  of  15.  The  last  column  provides  the  correlation  structure 
of  the  average  data.  The  estimated  correlations  remain  artifi¬ 
cially  high  in  this  case  because  averaging  reduces  the  varia¬ 
bility  of  the  data  about  the  year  cycle  which  intensifies  the 
artificial  increase  in  correlation  due  to  the  year  cycle.  The 
correlation  structure  revealed  in  Table  IV. B. 2  for  individual 
years  closely  resembles  that  of  an  AR(1)  model,  in  that  the 
k-step  correlation  is  approximately  the  one-step  correlation 
raised  to  the  kth  power.  The  correlations  in  the  table  have 
a  tendency  to  be  slightly  higher  than  the  theoretical,  calcu¬ 
lated  value,  but  the  agreement  is  reasonably  good  for  about 
ten  steps.  Beyond  that  point  the  correlations  are  kept  up  by 
the  year  cycle,  which  is  not  as  prominent  in  the  yearly  data 
as  it  is  in  the  averaged  data.  If  nothing  else  the  disparity 
between  the  two  correlations  is  evidence  of  the  existence  of 
a  trend  in  the  data. 

At  this  point  sufficient  information  is  available  to  de¬ 
termine  some  characteristics  of  the  general  form  of  the  model 
for  representing  the  wind  speed  data.  As  noted  above,  the 
correlation  structure  is  similar  to  that  of  an  AR(1)  process 
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rnu*  >r-»oo  >om  >r  <sj«n  hhooooo*®  >o^»rs>i>*r^r<*r»«o<30QO 

0(0(SiZ)r«r«r«ri*Nr»r<*r»Kr»sr*>-g<o<o^)'U^^^^4Jo<04>>o>0^ 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


ofNjmc^'Oo-ooooo'O'T  <ocof“o<7'-^o,‘oco>ouvnn^^comrsi'f  ^  i 

~«rufv<vfNro^'^'r'r>Tn'y>T‘f  j 

00000000300000000000000000000000  I 

•  ••••••••••••••'••••••••••••a***  I 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO  I 


fnco^in- koo^ on  co  in  m  f\JiAsf  o  *-<  ^-4-  tn^m 

r^a®  >o  -O  *•  -«irwM  ao  <0  ^  ®  eg-*  -•«■*  o^o  f-  ®m  m«>*t  >f  •>*  u> 

CD>flU\«fnmN(MH^H-4H^H-K3000000000000000 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


(DNinsf(<>NNr«i^HH^^^HHHOoOOOOOOOOOOO  I 

OOOOOOOOOOOOOOOOOOOOOOOOO I  I  I  r  I  I O  I 

p-'f-rvj^/Ni^ao^^or^P-.oiniAiAOO'Oi/Nvn't m-*Oco e-e-e-coO"**  I 

doooooooaoooooocoooooooooooooooo  i 

rvo*»o^O"^0*«,'Ocomr^a*co-o>OirMf\»fiAir\r-<x>c«Of*»W'tfO*^»Mro  l 
eo^in^rtrtNNN^H^coceoooooooHOoocoooo  i 

dododddododooddddododooddodoodoo  | 

?TjO'Oco<7*-^,co-«^t7‘'43fO*^®^-kri'frsj^OOO^O'-,^cgO-^^CMrMfO  I 

c0Mfl>trt^fMNH-i^HOOOOOOOOeOOOOOOOOOOO  | 
OOOOOOOOOOOOOOOOOO I O I OOOOOOOOOOO  | 
•r -^co'T*— *,0rNjcDr»inm^-*OOO-*— «f\Jrgf\jrMm(nmf\i*-«oo<7‘®  | 

oooaoooooooaoooooocooooooooooooo  i 

0>0't‘*#«Ar^-Aa3m-<cor'>oiAu'i>OU>  -t  >f'^'\jrnrn-t,<\irM— «-*.-4COO  | 
©flflN*^NH^HOOOOoOOOOOOOoOOOOOOOOO  | 

. . . .  I 

oooaoooooooooooooooocoaooooooo »  O  I 

fno°oosom<7*'Omoaoe«.>o«T  >r  mrsj-*ooo-*<\iror\j.-«oa'fw»nmin  i 

®r-tns*^fnr\jo4rgrg^**«*i-*— »-*-h«-«**OCOoO  j 

. . •••••••••••••••••••••••! 

OOOOOOOOOOOOOOOCOOOOOCOOOOOOOOOO  I 

mo  rgc®  cummin  | 

(XjMir^^^iNr^tvj^^^^H^^H^HOoooooooao  i 

odddddddoooododdoddodoocoooododo  I 

>^»r  <3  O' 'TO'O  uvn«T'^'f'T rorgrvr^^O -^000-^000  oj  | 
ffl^in^<nNrt-*-40«ooooooocooooooococooo  l 

ddcoddddodcddcddoddddo  i*o  i  *  i  do  i  do  | 

o^oiPu>'00-or\ja'^(^o^-'£nmArnc\j— | 

©•0©vj©nrg('JHHH^OOOCOOCOCCOOCOOOOOOO  I 

•.••••••’•••••••••••••••••.••••A  I 

OOOOOOOOCOOOOOOOOOOO I OOOOOOOOOOO  I 

mo®®  m*fmmmrgr,j*~‘>-<0'<7' tf'C'cr®  r»  m®  ®cro  i 

crr-tn  ^©©(N^r^-^^rHH^w^HHHCOOOOOCOOOOH  j 

coooooooeocoocoocoooeooocooooooo  I 

-*®f*»®OtnPgO,>*0<'JOO'(7*O®®«om^m<M<M*,«'"<OC7'O'®fte®^>»**  I 
®  Oin^'fmmNNNN«WHrtiM1MHHH^'-'-H.jOOOOCOO  1 

• . . . .  •  I 

ooeocooccooooocooocorooooooooooo  I 

'O  -OtfMA  IMO0'p»^'0'0^‘A 

B^\n.r,n("f'(rjH-4^N^'-<^-4-«HH^HH^ooOOCOO 

oooooococoooooocooccocoeooeoceoo 

eooocooccocooooccoccccoecocooooo 

(M®  \f\ -sUnN—*® r*«-  >nm  sj  UMfMn 't  r>^ lf\ 

«4ir^mNN-MHHHNH«»^ooccccoeccccccoo 

oooooooccoccoc occoccccoccor cocoo 


*t  m®?*®  if  C  m  ®  r»<o  a*  O  ro  «*  m  ®  f**  ®  C*  ©^«ni 

«-<  »m  ^ cv  j  (m  rg  rsi  «\j  r  j  <nj  r\i  oj  »vm^  m 


253 


with  positive  correlation  although  nuances  may  appear  as  the 
year  cycle  is  removed.  Hence,  an  AR(1)  model  should  be  used 
as  a  starting  point  for  the  construction  of  the  model  for 
wind  speed. 

In  addition,  the  cyclic  nature  of  the  process  must  be 
modeled.  This  can  be  done  with  either  an  additive  or  multi¬ 
plicative  model.  An  additive  model  might  have  a  structure  as 
follows. 


n 


=  u  +  e  , 
Mn  ln' 


(IV. B. 3) 


where  {XR}  is  the  time  series  under  consideration,  pn  is  a 

deterministic  function  of  n,  and  the  innovative  sequence 

{en}  is  a  stationary  sequence  of  random  variables.  In  the 

usual  model  this  stationarity  implies  that  the  marginal  vari- 
2 

ance  a  is  constant  and  the  correlations  only  depend  on  the 
lag  (i.e.,  p(Xn,Xn+k)  =  p(k)).  Using  the  same  definitions 
the  multiplicative  model  would  have  the  form 


X 


n 


un£n' 


(IV. B. 4) 


where  again  the  {en>  sequence  is  stationary  and  independent 
of  un*  A  characteristic  of  the  additive  model  is  that  the 
coefficient  of  variation  is  a  function  of  the  value  of  uR. 
The  multiplicative  model  produces  a  constant  coefficient  of 


variation.  Since  the  data  has  a  coefficient  of  variation  that 


is  essentially  constant,  in  the  crude  monthly  analysis,  the 
multiplicative  model  is  preferred. 

We  have  yet  to  determine  the  exact  form  of  the  mean,  u  , 
in  equation  IV. B. 4.  However,  we  do  know  that  this  mean  will 
have  a  yearly  cyclic  nature.  We  also  have  yet  to  determine 
the  general  structural  nature  of  the  innovative  process  e n  . 
These  subjects  are  addressed  in  the  following  sections. 

C.  THE  FORM  OF  THE  MEAN;  DETRENDING  THE  DATA 

Two  basic  models  were  considered  to  represent  the  mean. 

The  first  was  a  single  harmonic  sinusoidal  model 

un  =  a  +  bx  sinewy)  *  bj  cos(^%>  -  a  ♦  k  cos  +  8)  , 

(IV.C.l) 

2  2  1/2  -1  bl 

where  k  =  (b.^  +  b^)  and  e  =  tan  (-  g— )  .  The  second  was 
an  exponential  sine  with  one  harmonic 


n 


=  e 


.  i  •  f  2 Tin  %  .  i _ _ t  2un  « 

a+bisin (2$io) +b2cos  29T0 


=  e  e 


a  kcos(292O,0) 


(IV. C. 2) 


The  second  model  has  the  theoretical  advantage  that  it  can 
not  be  negative  and  will  represent  higher  harmonics  in  a  com¬ 
pact  form.  The  sinusoidal  model  may  or  may  not  be  negative 
depending  on  the  values  for  a,  b^,  and  b2 •  In  spite  of  the 
theoretical  preference  for  the  exponential  sin®,  both  models 
were  used  initially  to  see  if  either  produced  significantly 
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better  results.  Note  that  if  k  is  small  the  models  are 
hardly  distinguishable.  If  k  is  large,  the  exponential  sin 
will  clip  at  low  values  and  is  a  cycle  that  would  have  many 
harmonics  in  its  Fourier  transform . 

The  values  for  the  constants  in  equation  IV.C.l  were 
determined  in  a  straightforward  procedure  using  the  least- 
squares  regression  procedure  of  MINITAB  and  the  data  averaged 
over  15  years.  These  estimates  could  also  have  been  obtained 
from  the  periodogram  at  =  2ttN,  1(00^)  ~  {(b^)  c  +  (b^^), 
using  the  relations 

-  NX. 

a  -  l  =  x;  (IV. C. 3) 

i=l  N 

-  N  ,  . 

b.  =  2  l  X.  sin(^~ )/N  =  imaginary  component  (IV. C. 4) 

i=l  of  periodogram  at 

2tt/N  ; 

N  2  . 

b2  =  2  l  X .  cos  (~=-)  /N  =  real  component  of 

“  i=l  1  periodogram  at  2tt/N. 

2 

The  variance  of  these  estimates  is  2a  /N  if  the  X- 's  are 

e  1 

independent,  but  since  this  is  clearly  not  the  case  here, 
estimates  of  the  variance  of  the  estimates  cannot  be  obtained 
directly.  The  results  of  the  estimation  are  contained  in 
column  1  of  Table  IV.C.l. 

Similar  results  were  obtained  for  the  constants  in  IV. C. 2 
by  a  slightly  more  complicated  procedure.  In  order  to  use  a 
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TABLE  IV.C.l 


Parameter  Estimates  for  Models  of  the  Mean  Value 
Function  of  the  Wind  Velocity  Data 

ESTIMATE 

1  harmonic  1  harmonic  2  harmonic  2  harmonic  harmonic 

PARAMETER  sine  exp  sine  sine  exp  sine  exp  sine 

a/a'  10.230  2.309  10.230  2.307  2.307 

b±  -0.176  -0.011  -0.175  -0.011  -0.011 

b2  2.560  0.260  2.566  0.260  0.260 

b3  -  -  -0.593  -0.057  -0.057 

b.  -  -  -0.397  -0.054  -0.054 

4 

b5  -  -  -  -  0.014 

bg  -  -  -  -  0.001 

b?  -  -  -0.010 
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least  squares  approach,  a  linear  relationship  must  be  estab¬ 
lished  for  the  mean  value  of  the  process.  Taking  logs  is 
the  obvious  technique  to  employ,  but  this  introduces  a  compli¬ 
cation.  Taking  logs  and  expectation  of  IV. C. 2  we  have 


E  (In  X  )  =  In  u  +  E(ln  e  ) 

n  n  n 


*  *  +  »i  sin <$&)  +  b2  cosl^j)  +  c, 


For  example,  if  the  {en}  sequence  is  marginally  distributed 
as  a  unit  Gamma  variate,  G(l,k),  then  c  =  $  (k)  ,  where  >•  (k) 
is  the  digamma  function  (derivative  of  In  ;  (k) )  .  See  Cox  and 
Lewis  [Ref.  29],  pages  24-27.  The  value  of  the  constant  c 
will  be  combined  with  the  constant  a  in  the  least  squares 
estimation  using  the  In  Xn's,  giving  the  constant  a'  =  a+c. 

To  estimate  a+c  without  making  Gamma  assumptions  for  the 
innovative  process,  the  X^'s  are  divided  by 


=  bl  sin(^V  +  b2  cos(J^70) 


to  give  X^.  The  data  is  then  divided  by  the  average  of  the 

X„ ' s  which  estimate  e  ^a+c^ .  The  result  of  this  is  a  series 
n 

with  mean  value  (within  statistical  fluctuation)  of  1  if  the 
model  for  the  cycle  is  correct.  The  values  obtained  are  listed 
in  column  2  of  TAble  IV.C.l.  The  results  of  these  estimates 
are  in  Figure  IV.C.l.  In  this  figure  the  average  data 


258 


is  plotted  against  the  value  computed  for  yn  using  both  of 
the  models  under  consideration. 

When  using  a  multiplicative  model,  the  residuals  are 
formed  by  dividing  the  raw  data  by  the  mean.  The  results  of 
this  procedure  are  presented  in  Figures  IV.C.2a  through  IV.C.2p 
using  the  exponential  sine  model  for  the  mean.  The  results 
are  not  significantly  different  using  the  sinusoidal  model 
of  the  means.  Hence,  only  the  results  for  the  average  data 
are  presented  for  this  case  in  Figure  IV. C. 3. 

The  log  periodogram  of  the  average  data  detrended  using 
the  sinusoidal  model  for  the  mean  is  shown  in  Figure  IV. C. 4. 

A  five-step  moving  average  of  this  log  periodogram  is  pre¬ 
sented  in  Figure  IV. C. 5.  The  detrending  has  clearly  reduced 
the  importance  of  the  yearly  cycle,  but  still  shows  some  evi¬ 
dence  of  a  six  month  cycle  and  six  and  twelve  hour  cycles. 
Similar  information  is  provided  for  the  average  data  detrended 
using  the  exponential  sine  model  for  the  mean  in  Figures  IV.C.6 
and  IV. C. 7.  This  model  does  not  reduce  the  effect  of  the 
yearly  cycle  as  much  as  the  sinusoidal  model  for  the  mean, 
but  still  shows  the  six  month  cycle  as  being  important  and 
some  evidence  of  six  and  twelve  hour  cycles. 

Since  the  exponential  sine  has  the  theoretical  advantage 
of  being  non-negative  and  both  models  of  the  mean  produce 
similar  results  when  applied  to  the  data,  the  exponential 
sine  is  selected  as  the  model  of  choice  and  the  analysis  is 
continued  using  it  exclusively. 
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*  >«)ur<‘  IV.  <’.6.  fH’i  shows  dominance*  of  six-month  cycle  and  |)H*s«*n<*o  of  six  and  twelve  hour  cycles  aft 

t  (-niov.il  of  yearly  cycle. 
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D.  RESIDUAL  PROCESS  PROBABILITY  STRUCTURE 

Having  removed  the  dominant  seasonal  effect  from  the 
data,  it  is  possible  to  investigate  the  structure  of  the 
residual  process  t  in  the  model 


X 


n 


ynen- 


(IV.D.l) 


The  two  facets  of  the  probability  structure  of  the  stationary 
process  {en}  which  were  addressed  in  Chapter  II  were  the 
marginal  distribution  and  the  correlation  structure.  The 
residuals  produced  by  dividing  the  raw  data  by  the  appropriate 
value  of  the  mean  were  supplied  to  HISTF ,  a  histogram  and  box 
plot  routine  developed  at  the  Naval  Postgraduate  School. 
Histograms  for  each  year  and  the  entire  data  set  were  pro¬ 
duced.  These  histograms  are  presented  in  Figures  IV.D.la 
through  IV.D.lp.  The  shape  of  the  histograms  is  consistent 
over  the  years  and  indicates  that  a  Gamma  distribution  is 
appropriate  for  modeling  the  innovative  factors.  The  param¬ 
eter  k  can  be  estimated  as  the  reciprocal  of  the  coefficient 
of  variation  squared  (see  equation  II. B. 4. 16).  The  estimated 
value  of  k  for  each  year  is  given  in  Table  IV.D.l. 

A  careful  examination  of  the  statistics  associated  with 
the  histogram  will  reveal  that  the  values  for  the  skewness 
and  kurtosis  are  low  compared  to  the  theoretical  values  for 
the  Gamma  distribution,  namely  2//k  and  6/k,  respectively. 
However,  this  is  not  unexpected  when  one  recalls  the 
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Histoyraa  and  boxplot  of  1957  detrended  data 
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Figure  XV.D.le.  Histogram  and  boxplot  of  1959  detrended  date. 
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Figure  IV.D.lh.  Histogram  and  boxplot  of  1962  detrended  data 
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figure  IV. D. In.  Ulitogrea  end  boxplot  of  196?  detrended  date 


Table  IV.D.l 


Moment  Estimate  of  the  Gamma  Shape  Parameter  by  Year 


Year 

Estimate 

1955 

3.96 

1956 

4.16 

1957 

4.38 

1958 

3.68 

1959 

3.65 

1960 

4.45 

1961 

3.87 

1962 

3.87 

1963 

3.86 

1964 

4.49 

1965 

4.32 

1966 

4.04 

1967 

5.08 

1968 

4.24 

1969 

5.41 
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discretization  and  clipping  that  has  apparently  occurred 
while  the  data  was  collected  and  processed.  As  a  rough 
check  on  the  extent  of  the  clipping  required  to  produce 
values  for  the  kurtosis  and  skewness  similar  to  those  for  the 
data  the  following  procedure  was  examined.  A  sample  of  2920 
values  for  a  Gamma  distribution  were  produced  with  mean  1  and 
shape  parameter  k  =  4.0  by  a  call  to  the  NPS  random  number 
generator  LLRANDOMII  (SGA11A)  .  A  histogram  of  these  values  was 
produced  sequentially  for  the  following  cases.  The  data  was 
clipped  so  that  all  values  over  four  were  set  equal  to  four, 
all  values  over  three  were  set  equal  to  three,  and  finally  the 
highest  ten  percent  of  the  data  was  set  equal  to  the  value  of 
the  2890  sample  order  statistic.  The  first  four  central  moments 
were  estimated  under  each  of  the  conditions.  The  results  are 
presented  in  Table  IV.D.2.  The  results  indicate  that  a 
clipping  of  the  top  ten  percent  of  the  data  will  yield  results 
for  skewness  and  kurtosis  comparable  with  those  observed  in 
the  data. 

TABLE  IV.D.2 


Mean 

SD 

CV 

Skewness 

Kurtosis 

Gamma 

0.986 

0.504 

0.511 

1.128 

2.114 

Cut  at  4.0 

0.986 

0.504 

0.511 

1.128 

2.114 

Cut  at  3.0 

0.984 

0.498 

0.506 

0.994 

1.165 

10%  Cut 

0.982 

0.489 

0.498 

0.861 

0.516 
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The  conclusion  that  the  innovative  factors  can  be  modeled 
as  random  variables  with  Gamma  marginals  when  combined  with 
the  conclusions  from  IV. B  specify  the  general  form  of  the 
model  to  be  possibly  that  of  the  GLAR ( 1 )  process,  although 
further  detrending  might  indicate  that  the  more  general 
GLARMA (p, q)  model  of  Chapter  II  might  have  to  be  used. 

Since  the  estimated  correlations  o (k)  are  affected  by 
remaining  trend  (as  seen  in  Table  IV.0.3-),  it  is  best  to 
examine  the  structure  of  the  dependency  process  via  the 
periodogram. 

Figures  IV. D. 2  through  IV. D. 5  show  the  periodogram  and 
log  periodogram  for  the  1955  and  1969  data  detrended  by  the 
single,  yearly  harmonic  exponential  sine  (see  equation  IV. C. 2) 
Superimposed  over  these  plots  is  the  spectral  density  and  log 
spectral  density  of  a  theoretical  AR(1)  process  with  correla¬ 
tion  equal  to  0.849  (see  Table  IV.D.3).  This  spectral  density 
is  also  the  spectral  density  of  the  GLAR ( 1 )  process.  We  have 

2 

f(w)  =  — - £-J -  ,  0  <  u  <  it,  (IV. D.  2) 

1  +  p  -  2  cos (w) 

with  o  -  0.849. 

All  of  these  plots  show  that  the  detrending  has  reduced 
the  importance  of  the  yearly  cycle  and  that  a  six  month  cycle 
has  now  become  the  dominant  factor.  The  theoretical  GLAR(l) 
spectral  density  fits  well  for  the  periodogram  after  the 
point  representing  the  six  month  cycle.  The  six  month  cycle 
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DATA:  3 -HOURLY  MIND  VELOCITIES;  OETRENO INGS  I  HARMONIC  EXPONENTIAL  SINE 
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process  with  correlation  of  O.M<>  is  super  imposed  over  pel  ioihujrsm  of  1  harmonic 


V  i«jurc*  I V .  O .  i  .  Sjjuc 

Hot  r 


has  now  become  the  dominant  factor.  Lacking  in  these  plots  is 
an  indication  of  a  time  of  day  effect.  The  appearance  of  a 
time  of  day  effect  is  limited  to  those  plots  which  use  the 
average  data  (see  IV.C.4  through  IV.C.7).  Figures  IV.C.4 
through  IV.C.7  and  IV. D. 2  thorugh  IV. D. 7  indicate  that  a 
further  refinement  in  the  model  of  the  mean  to  include  a  six 
month  cycle  may  be  helpful.  This  topic  is  covered  in  the 
next  section. 

The  correlation  structure  of  the  detrended  data  is  depicted 
in  Table  IV. D. 3  and  Figures  IV.D.6a  through  IV.D.6p  and  IV. D. 7. 
Since  the  one-harmonic  year  cycle  in  the  data  has  been  reduced, 
the  correlation  of  the  average  data  in  Table  IV. D. 3  more  closely 
reflects  that  of  the  average  of  the  fifteen  yearly  correlograms . 
The  higher  values  for  the  correlations  in  the  average  data  and 
its  failure  to  fall  below  0,20  may  be  indications  that  a  trend 
still  exists  in  the  data  (the  six  month  trend)  which  is  arti¬ 
ficially  inflating  these  values.  This  may  be  a  further  indi¬ 
cation  of  the  desirability  of  including  further  cycles  in 
the  model  of  the  mean.  The  slight  increases  in  the  correlations 
for  lags  of  8,  16,  and  24  in  the  average  correlogram,  Figure 
IV. D . 6p,  and  the  correlogram  for  the  average  data,  Figure 
IV. D. 7,  may  indicate  a  small  time  of  day  effect. 

E.  REFINING  THE  FORM  OF  THE  MEAN;  A  FURTHER  DETRENDING 

Since  several  plots  in  the  previous  section  indicated  that  a 
six  month  cycle  had  become  the  dominant  factor  since  the  removal 
of  the  one-harmonic  year  cycle,  a  further  refinement  for  the  model 


307 


CORRELATION  LAG 


.CORRELATION  LAG 


CORRELATION  LAC 


CORRELATION  LRG 


CORRELATION  LAG 


.  CORRELATION  LAG 


# 


CORRELATION  LAG 


•CORRELATION  LRG 


•CORRELATION  LAG 


CORRELATION  LAG 


of  the  mean  was  considered.  In  this  refinement,  the  sinu¬ 
soidal  model  for  the  mean  was  reintroduced  to  the  analysis 
to  see  if  the  addition  of  the  second  cycle  would  allow  this 
model  to  outperform  the  exponential  sine .  With  the  two 
cycles  included  the  sinusoidal  model  becomes 

Un  =  a  +  b1  sm  ( 2gTJo’)  +  b2  (2920‘)  +  b3  Sin  (TT^) 

+  b4  cos  -  (IV.E.l) 

The  exponential  sine  becomes 


=  e 


a+b^sin ( 


2  7T  n  >  { .  ,  2  tt  n  >  ,  *  >  , 

2920  +b2  2920  +fa3sin ( 


27rn 

1460 


) +b .cos  ( 


2irn 

1460 


(IV. E. 2) 


Parameters  for  these  models  were  determined  following  the 
procedures  in  IV. C  and  the  estimated  values  are  listed  in 
Table  IV.C.l.  The  plots  of  the  two  resulting  values  for  the 
mean  are  presented  in  Figure  IV.E.l.  Since  the  two  curves 
are  nearly  identical  and  the  exponential  sine  is  preferred 
on  a  theoretical  basis,  the  sinusoidal  model  was  not  further 
considered. 

Figures  IV.E.2  through  IV.E.5  show  the  periodogram  and 
log  periodogram  for  1955  and  1969  after  detrending  with  the 
two  harmonic  exponential  sine.  The  value  for  an  AR(1)  proc¬ 
ess  is  superimposed  as  before  with  a  correlation  of  0.794 
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comparable . 


(see  equation  IV. D. 2) .  All  of  these  plots  show  the  yearly 
and  six  month  cycles  much  reduced  in  importance.  They  also 
show  some  weak  time-of-day  effects,  but  these  are  more 
noticeable  in  the  1955  data. 

The  correlation  structure  of  the  data  is  shown  in  Table 
IV.E.l.  The  average  data  correlations  are  now  lower  than 
those  of  the  average  of  the  fifteen  yearly  correlations. 

They  also  drop  more  quickly  than  that  of  the  average  of  the 
fifteen  yearly  correlations  and  eventually  go  below  zero. 

This  is  another  indication  that  the  trend  in  the  average  data 
has  been  largely  removed.  Figures  IV. E. 6  and  IV. E. 7  show  the 
periodogram  and  log  periodogram  for  the  averaged  data.  As 
has  been  noted  previously,  the  time  of  day  effects  are  more 
noticeable  in  the  averaged  data  than  they  are  in  the  data 
for  a  single  year.  However,  the  effects  are  prominent  enough 
to  warrent  further  consideration.  This  subject  will  be  ad¬ 
dressed  in  Section  IV. G. 

F.  RESIDUAL  ANALYSIS 

Since  a  first-order  autoregressive  model  appears  to  be  a 
good  fit  to  the  innovation  process  {sn},  we  need  to  examine 
this  hypothesis  more  clearly.  If  we  were  dealing  with  a 
linear  AR(1)  model  for  the  residual  process 

(IV.F.l) 

n  n-l  n 

where  Yn  is  a  sequence  of  iid  random  variables,  then  computing 
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Prominent  time  of  day  cycles  are  evident  in  the  averaqe  data  after  2  hat  monte  detrendlnq. 


(IV. F. 2) 


as  an  estimate  of  Y„  would  be  of  interest.  The  estimated 

n 

(YnJ  should  have  a  flat  spectrum.  Using  the  Gamma  generation 

scheme  of  equation  II. A. 5  (i.e.,  e  =  B  z  ,  +  C  G  )  reduces 

n  n  n-1  n  n 

the  value  of  differencing  since  the  coefficients  in  the  genera¬ 
tion  scheme,  Bn  and  Cn,  are  continuous  random  variables  and 
not  constants.  However,  this  differencing  procedure  may  pro¬ 
duce  some  insight  to  the  data.  Hence,  the  differences 


A  A 

Yn  =  en  "  p(1)cn-l  (IV. F. 3) 

were  produced,  where  p  is  a  one-step  (lag  one)  correlation 

for  the  two  harmonic  detrended  data  and  e  and  e  .  are  two 

n  n-l 

harmonic  detrended  data  values. 

Since  the  data  has  been  detrended  and  the  differencing 
serves  to  remove  the  dependence  from  the  data,  one  would  expect 
the  periodogram  of  the  detrended,  difference  data  to  be  flat. 

The  periodogram  and  log  periodogram  for  the  detrended,  dif¬ 
ferenced  data  are  presented  in  Figures  IV.F.l  and  IV. F. 2. 

With  the  exception  of  a  relatively  strong  indication  of  a 
six  and  twelve  hour  cycle,  the  periodogram  is,  in  fact, 
reasonably  flat.  The  log  periodogram  indicates  the  same 
characteristics . 

The  correlogram  for  the  detrended,  differenced  data  is  Figure 
IV. F. 3.  There  are  two  key  points.  First,  all  the  correlations 
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iMt.i  |<ic*whi to  remove  the  autorctjresfl ivo  component  shows  presence  of  six  and  twelve  hour  cycles 
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are  relatively  low,  indicating  that  the  dependence  structure 
has  been  largely  removed.  Second,  the  alternation  of  the 
sign  of  the  correlations  is  an  indication  that  there  still 
exists  an  important  cyclic  component  in  the  data,  in  particu¬ 
lar  an  alternation  of  twelve,  or  six  hours.  Differencing 
two  sine  functions  (i.e.,  sin(^|— )  -  s in ( —  - — )  will  pro¬ 
duce  a  cycle  with  period  of  two  if  they  have  non-zero  ampli¬ 
tude.  Therefore,  the  alternation  of  the  correlations  is 
evidence  that  an  important  cycle  still  remains  in  the  data. 

G.  A  FURTHER  REFINEMENT  OF  THE  MEAN:  THE  LAST  DETRENDING 

Since  the  evidence  of  the  preceding  two  sections  suggests 
that  there  is  still  one  important  cycle  in  the  data,  a  further 
refinement  of  the  model  for  the  mean  was  undertaken.  The 
evidence  suggests  that  there  may  be  six  and  twelve  hour  cycles 
in  the  data.  These  cycles  may  be  the  result  of  the  passage 
of  pressure  fronts  over  the  data  collection  location. 

Only  the  exponential  sine  model  for  the  mean  was  con¬ 
sidered  in  the  final  detrending.  The  final  model  for  the  mean 
was 

yn  =  EXP [ a+b1 sin ( 2 ^ )  + b2cos  ( jqJq)  +  b2sin{1460*  +  b4cos (14g0) 
+  b4sin(^£)  +  bgCOs  +  b?cos  <-—)  ]  (IV.G.l) 

One  should  note  that  the  sine  function  with  a  period  of  two 
is  omitted  from  the  model.  This  is  because  the  sin  (mr)  is 
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identically  zero  if  n  is  an  integer.  The  implication  of  this 
is  that  we  essentially  lose  the  ability  to  determine  the 
phase  shift  for  the  cycle  with  period  two.  This  may  mean 
that  our  attempt  to  remove  the  six  hour  cycle  will  not  be 
completely  successful.  Parameters  for  the  model  in  IV.G.l 
were  produced  by  the  same  techniques  used  previously  (see 
Section  IV. C  and  Table  IV.C.l) . 

Figures  IV.G.l  and  IV.G.2  show  the  periodogram  and  log 
periodogram  for  1955  data  after  detrending  using  the  model 
of  the  mean  in  IV.G.l  (see  also  Table  IV.C.l).  With  the 
exception  of  the  six  hour  cycle,  the  periodogram  compares 
favorably  with  the  theoretical  AR(1)  periodogram  superimposed 
over  it.  The  log  periodogram  shows  the  same  characteristics. 
Similar  information  is  presented  for  1969  in  Figures  IV. G. 3 
and  IV. G. 4.  The  strength  of  the  six  hour  cycle  is  reduced 
for  this  year.  Finally,  the  periodogram  and  log  periodogram 
for  the  averaged  data  are  presented  in  Figures  IV. G. 5  and 
IV. G. 6.  The  comparison  of  the  averaged  data  with  a  theoreti¬ 
cal  AR ( 1 )  process  is  considered  acceptable. 

Note  too  that  in  Table  IV.G.l  the  15  year  average 
correlation  is  commensurate  with  the  correlation  computed 
from  the  averaged  data.  Thus  the  discrepancy  between  these 
quantities  noted  in  Table  IV. B. 2  has  been  explained. 

It  may  be  worth  noting  in  passing  that  a  surprising  result 
of  this  analysis  is  the  failure  to  detect  any  multiple-day 
cycles.  Apparently  some  meteorologists  believe  that  there 
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hour  cycle  after  4  harmonic  dot  rend in<j  imliratos  some  hotor  cn|f,no  i  t  y  among  yea*  a. 


iijure  IV 


Figure  iv. 5.  Average  data  shows  autoregressive  behavior  with  no  cycles  after  ♦  harmonic  detrending. 
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riquro  IV. 


is  a  six-day  weather-cycle  driven  by  the  passage  of  storms. 
This  analysis  has  failed  to  detect  any  such  cycle.  It  may 
be  that  the  high  correlation  among  the  data  and  the  expecta¬ 
tion  that  the  actual  storm  cycle  will  be  reflected  in  the 
data  has  created  an  impression  that  these  cycles  exist  in 
the  data  when,  in  fact,  they  do  not.  This  confusion  of 
quasi-cycles  produced  by  high  positive  correlation  and  com¬ 
pletely  deterministic  cycles  is  common  in  applied  science. 
Figure  IV. G. 7  shows  a  sample  path  for  a  GLAR(l)  process  with 
high  correlation,  p (1)  =  0.85.  Although  it  may  be  tempting 
to  conclude  that  this  process  is  showing  evidence  of  a 
cyclic  nature,  there  is  no  deterministic  cycle  in  the  data 
shown.  The  behavior  displayed  in  this  figure  is  typical  of 
an  autoregressive  process  with  high  correlation,  and  no  cycle. 

A  table  of  correlations  for  the  4  harmonic  detrended 
data  is  provided  in  Table  IV.G.l.  Its  characteristics  are 
much  the  same  as  those  of  the  two  harmonic  detrended  data. 

H.  SUMMARY 

The  model  suggested  for  the  representation  of  the  wind 
speed  data  now  has  the  following  form.  The  basic  structure 
is  that  of  a  multiplicative  model,  that  is  it  has  the 
form 

xn  =  »nc-n>  n  ■  i/2,...  .  (IV.H.l) 
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The  {Xn>  sequence  represents  the  raw  wind  speed  data.  The 
Un  is  a  deterministic  function  of  n.  The  innovative  terms 
{ }  are  modeled  by  a  GLAR(l)  process. 

The  GLAR(l)  process  was  discussed  and  analyzed  in  Section 
II. B.  The  generation  scheme  presented  in  equation  II. B. 1.1 
is  repeated  here  (with  ei  replacing  . 


n 


=  A  e  ,  + 
n  n-l 


B  G  . 
n  n 


(IV. H. 2) 


The  innovative  sequence  {en>  is  itself  correlated.  The 
parameters  of  the  GLAR(l)  process  control  the  correlation 
structure  of  the  model.  (See  Section  II. B. 2,  in  particular 
equation  II. B. 2.1.) 

The  mean  un  has  been  modeled  as  a  four  harmonic  exponen¬ 
tial  sine  function. 


“n  ‘  Exp  (a+b^in  <$&>  +b2cos  (^)  +b3sin  <$&>  *b4co,  <$&> 


+b5sin  (2j£)  +b6cos  (2j£)  +b7cos  (^)  ] 


(IV. H. 3) 


The  four  harmonics  represented  are  a  yearly  cycle  (coeffi¬ 
cients  b^  and  b2)  >  a  six  month  cycle  (coefficients  b3  and 
b4) ,  a  twelve  hour  cycle  (coefficients  bg  and  bg)  and  a  six 
hour  cycle  (coefficient  b^) .  The  values  for  these  parameters 
and  the  "a"  parameter  are  found  in  Table  IV.C.l. 

The  innovative  terms  are  modeled  by  the  GLAR(l)  process. 
The  parameter  values  for  k  and  q  were  determined  to  be  2.843 
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and  0.727,  respectively,  by  using  the  numerical  approximation 
to  the  maximum  likelihood  method  described  in  Section  II. B. 4. 
The  data  used  for  this  evaluation  were  the  residuals  produced 
by  the  single  harmonic  exponential  sine  model  of  the  mean. 

The  parameters  were  not  recomputed  for  the  two  or  four  har¬ 
monic  exponential  sine  model  because  of  time  limitations. 

These  parameter  values  give  a  correlation  of  0.744.  This 
is  somewhat  less  than  the  average  correlation  of  0.826  for 
the  single  harmonic  residuals  (see  Table  IV.E.l) .  However, 
this  deviation  is  not  considered  serious.  This  is  because 
the  estimates  produced  by  the  four  harmonic  detrended  data 
may  differ  from  those  produced  by  the  two  harmonic  data  and 
the  correlations  for  the  one  harmonic  data  are  modified  by 
the  presence  of  the  six  month,  twelve  hour  and  six  hour  cycles. 

The  simulation  study  of  Section  II. E  indicated  that  for 
large  k  and  high  correlation  the  standard  deviation  of  the 
maximum  likelihood  estimates  was  about  half  that  for  the 
moment  estimates  (see  Figure  II.E.l  and  Table  II. E. 2).  In 
addition,  neither  estimation  procedure  had  any  apparent  bias. 
For  these  reasons  the  maximum  likelihood  estimates  are  pre¬ 
ferred  over  the  moment  estimates  in  this  case  unless  com¬ 
puter  time  is  limited. 
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